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ABSTRACT 

We present physical results obtained from simulations using 2+1 flavors of domain wall quarks 
and the Iwasaki gauge action at two values of the lattice spacing a, (a~ l = 1.73 (3) GeV and 
<3 _1 = 2.28 (3) GeV). On the coarser lattice, with 24 3 x 64 x 16 points (where the 16 corresponds 
to L s , the extent of the 5 th dimension inherent in the domain wall fermion (DWF) formulation 
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of QCD), the analysis of ref. OJ] is extended to approximately twice the number of configura- 
tions. The ensembles on the finer 32 3 x 64 x 16 lattice are new. We explain in detail how we use 
lattice data obtained at several values of the lattice spacing and for a range of quark masses in 
combined continuum-chiral fits in order to obtain results in the continuum limit and at physical 
quark masses. We implement this procedure for our data at two lattice spacings and with uni- 
tary pion masses in the approximate range 290-420 MeV (225-420 MeV for partially quenched 
pions). We use the masses of the % and K mesons and the Q. baryon to determine the physical 
quark masses and the values of the lattice spacing. While our data in the mass ranges above are 
consistent with the predictions of next-to-leading order SU(2) chiral perturbation theory, they are 
also consistent with a simple analytic ansatz leading to an inherent uncertainty in how best to per- 
form the chiral extrapolation that we are reluctant to reduce with model-dependent assumptions 
about higher order corrections. In some cases, particularly for f n , the pion leptonic decay con- 
stant, the uncertainty in the chiral extrapolation dominates the systematic error. Our main results 
include f n = 124(2) sta t(5) S yst MeV, fa/ fx = 1.204(7)(25) where fx is the kaon decay constant, 



m 



(2GeV) = (96.2 ±2.7) MeV and m^f (2 GeV) = (3.59 ±0.21) MeV (m s /m ud = 26.8 ± 1.4) 
where m s and m uc i are the mass of the strange-quark and the average of the up and down quark 
masses respectively, [E MS (2GeV)] 1 / 3 = 256(6) MeV, where E is the chiral condensate, the Som- 
mer scale r = 0.487(9) fm and r x = 0.333(9) fm. 
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I. INTRODUCTION 



For several years now, the RBC and UKQCD Collaborations have been undertaking a major pro- 
gramme of research in particle physics using lattice QCD with Domain Wall Fermions (DWF) 
and the Iwasaki gauge action. In the series of papers til-yD, we studied general properties of en- 
sembles with an inverse lattice spacing of a -1 = 1.73(3) GeV (corresponding to /3 = 2.13) and 
with unitary pion masses m n > 330 MeV (partially quenched m n > 240 MeV). The number of 
points in these ensembles are 16 3 x 32 x 8 M\, 16 3 x 32 x 16 |Q] and 24 3 x 64 x 16 Jl], where the 
fifth dimension is a feature of DWF and is not visible to low-energy physics which remains four- 
dimensional. We do not review the properties of DWF here, beyond underlining their physical 
chiral and flavor properties which we exploit in much of our wider scientific programme. We have 
used these ensembles to investigate a broad range of physics, including studies of the hadronic 
spectrum, mesonic decay constants and light-quark masses yj], the evaluation of the Bv parameter 
of neutral-kaon mixing |1, 4], the calculation of the form-factors of decays [Q, studies in 
nucleon structure O-la] and proton decay matrix elements 
study of the masses and mixing of the r\ and r\' mesons 



,10] and very recently the first lattice 
Hill as well as a determination of the 



matrix elements relevant for neutral 5-meson mixing in the static limit J12Q. A key limiting factor 
in the precision of these results was that the simulations were performed at a single lattice spacing. 
In this paper we remove this limitation, by presenting results for the spectrum, decay constants and 
quark masses obtained with the same lattice action using ensembles generated on a 32 3 x 64 x 16 
lattice at a second value of the lattice spacing corresponding to /3 =2.25, for which we will see 
below that cT x = 2.28(3) GeV. Now that we have results for the same physical quantities with the 
same action at two values of the lattice spacing we are able to perform a continuum extrapolation 
and below we will present physical results in the continuum limit. 

Since the most precise results at /3 = 2.13 were obtained on the 24 3 x 64 x 16 H 1 fj lattices, as a 
shorthand throughout this paper we will refer to these lattices as the 24 3 ensembles and label the 
new lattices at /3 = 2.25 as the 32 3 ensembles. 

The new 32 3 ensembles at /3 = 2.25 will, of course, be widely used also in our studies of other 
physical quantities. In this first paper however, we discuss their properties in some detail (see 
Sec. UIl). In this section we also discuss reweighting which allows us to eliminate one source of 
systematic uncertainty. While at present we cannot simulate with physical u and d quark masses, 
there is no reason, in principle, why we cannot simulate with the physical strange quark mass. 



4 



The difficulty however, is that we don't know a priori what this mass is and so in practice the 
simulations are performed with a strange quark mass which is a little different from the physical 
one. As explained in Section llTDl the technique of reweighting allows us to correct a posteriori 
for the small difference in the simulated and physical strange quark masses. In SectionUm we 
present updated raw results for the pion and kaon masses and decay constants and the mass of 
the ^-baryon on the 24 3 ensembles which have been extended beyond those discussed in ref.yj]. 
SectionHVl contains the corresponding results on the 32 3 ensembles. In these two sections we also 
present the raw results for the masses of the nucleon and A baryons from the two ensembles, but in 
contrast to the mesonic quantities a description of their chiral behaviour and extrapolation to the 
continuum limit are postponed to a future paper. 

The price we pay for using a formulation with good chiral and flavor properties is the presence 
of the fifth dimension and the corresponding increase in computational cost. The lightest unitary 
pion which we have been able to afford to simulate has a mass of 290 MeV and so, in addition 
to the continuum extrapolation we need to perform the chiral extrapolation in the quark masses. 
In Sec.|V]\ve present a detailed explanation of how we combine the chiral and continuum extrap- 
olations in an attempt to optimize the precision of the results, exploiting the Symanzik effective 
theory approach as well as chiral perturbation theory and other ansatze for the mass dependence of 
physical quantities. Having explained the procedure, we then proceed in Section lVEl to discuss the 
results, to determine the physical bare masses and lattice spacings as well as to make predictions 
for the pion and kaon decay constants. In particular we find that the ratio of kaon and pion decay 



constants |78|] 



fj- = 1.204 ±0.026, (1) 
fn 

where the error is largely due to the uncertainty in the chiral behaviour of f n as explained in 
Sec. lVE3l From the chiral behaviour of the masses and decay constants we determine the corre- 
sponding Low Energy Constants (LECs) of SU(2) Chiral Perturbation Theory (ChPT). 
Among the most important results of this paper are those for the average u and d quark mass and 
for the strange quark mass which are obtained in SecJVTl 

rojjjf (2GeV) = (3.59 ±0.21) MeV and mp(2GeV) = (96.2 ±2.7) MeV. (2) 

The masses are presented in the MS scheme at a renormalization scale of 2 GeV, after the renor- 



malization to symmetric momentum schemes has been performed non-perturbatively 



the conversion to the MS scheme has been done using very recent two-loop results [15 



13 



m. 



)M and 



5 



Section lVIIl contains a discussion of the topological charge and susceptibility of both the 24 3 and 
32 3 ensembles and in Sec lVIIII we summarise our main results and present our conclusions. There 
are three appendices. Appendix [A] contains the chiral extrapolations performed separately on 
the 24 3 and 32 3 ensembles. This is in contrast with the procedure described in Section lVEl in 
which the chiral and continuum extrapolations were performed simultaneously with common fit 
parameters at the two spacings. Appendix |B] contains a detailed analysis of a subtle issue, the 
normalization of the partially conserved axial current. For domain wall fermions this is expected 
to deviate from the conventionally normalized continuum current by terms of order am ies , where a 
is the lattice spacing and m res is the residual mass Jl, 17]. Current simulations are now becoming 
sufficiently precise that these effects need to be understood and quantified and the method proposed 
in appendix|B] in which the 0(am res ) effects are absent, is implemented in the numerical analyses 
throughout the paper. Finally Appendix O contains a discussion of the expected statistical errors 
when reweighting is performed on Monte Carlo data to obtain results with a different action from 
that used to generate the data. 

We end the Introduction with an explanation of our notation for quark masses [1]. When discussing 
unitary computations, with the valence and sea quarks degenerate, we call the bare light (u or d) 
quark mass m/ and the bare heavy (strange) quark mass m^. m uc j and m s refer to the physical values 
of these masses (we work in the isospin limit so that the up and down quarks are degenerate). For 
the partially quenched computations we retain the notation m/ and m/, for the sea-quark masses, 
but use m x and m y for the valence quarks. A tilde over the mass indicates that the residual mass 
has been added, m q = m q + m tes ; it is m which is multiplicatively renormalizable. 



II. SIMULATION DETAILS AND ENSEMBLE PROPERTIES 



As described in Ref. 



. HE, 



180, we generate ensembles using a combination of the DWF formula- 



tion of Shamir IU9I1 and the Iwasaki gauge action 112011 . For the fermionic action we use a value of 
1.8 for the "domain wall height" M5 and an extension of the 5 th dimension of L s = 16. In addition 
to the new ensembles generated on a 32 3 x 64 lattice volume and a gauge coupling /3 = 2.25, 
we have also significantly extended the 24 3 x 64, /3 =2.13 ensembles generated in our previous 
study [1]. As indicated in Tab. Uwe have extended the m/ = 0.005, 24 3 x 64 ensemble from 4460 
to 8980 MD units while the m/ = 0.01 ensemble has been extended from 5020 to 8540 MD units. 
The three 32 3 x 64 ensembles that are first reported here are also shown in Tab. IJ and those with 
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light quark masses of 0.004, 0.006 and 0.008 contain 6856, 7650 and 5930 MD units respectively. 



A. Ensemble Generation 



For the generation of both the 24 3 x 64 and 32 3 x 64 ensembles, we employ the "RHMC II" 
algorithm described in Ref. [1]. More specifically, the simulation of two light quarks and one 
strange quark is carried out using a product of three separate strange quark determinants each 
evaluated using the rational approximation. The 2 flavors of light quarks are preconditioned by 
the strange quark determinant j2l\ \. While the preconditioning mass does not have to be the same 
as the strange-quark mass, we found that the strange-quark mass is close to being optimal in DWF 
simulations in tests on smaller volumes. 

Using the notation 0(m/) = D^ DWF (Ms,mi)DowF(M^,mi), the fermion determinant including the 
contribution from the Pauli-Villars fields and evaluated on a fixed gauge configuration can be 
written as 



det 



0(1)3/2 



det 
det 



0(m,) 



3/2 



det 



3>{m s ) 



det 



0(1) 



det 



0(m,; 

0(1) 



det 



0(m,) 



(3) 
(4) 



In the third line we explicitly show how this ratio of determinants is implemented using the ratio- 
nal approximation. Here M a {x) denotes x° evaluated using the rational approximation and each 
determinant is evaluated using a separate set of pseudofermion fields. An Omelyan integrator Il22n 
with the Omelyan parameter A = 0.22 was used in each part of evolution. 

Given the disparate contributions to the molecular dynamics force coming from the gauge action 
and the different factors in Eq. © we follow the strategy of Ref. Il23n and increase performance by 
simulating these different contributions with different molecular dynamics time step granularities. 
In particular, the suppression of the force from the light quark determinant that results from the 
Hasenbusch preconditioning allows us to evaluate the computationally expensive force from the 
light quark using the largest time step among the different terms, decreasing the computational cost 
significantly. As a result, we divide our simulation in such a way that Alight : A^eavy : A? gau g e = 1 : 
1:1/6 which gave a good performance, measured in flops per accepted trajectory in tuning runs 
performed separately. (Note, the nature of the Omelyan integrator makes htheavy effectively half 
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of Atughf.) This ratio of time steps was used for all the ensembles studied here. However At[ ig i lt 
was varied from ensemble to ensemble to reach an approximate acceptance of 70%. The precise 
numbers that were used are listed in Tab. HI 

In addition, we chose to simulate with a trajectory length T = 2 for the 32 3 ensembles, twice 
that used for the 24 3 ensembles. While a longer trajectory length may be expected to reduce the 
autocorrelation between configurations, the time for a trajectory scales very nearly linearly in the 
trajectory length. In comparisons between x = 1 and x = 2 trajectory lengths we were not able 
to recognize any statistically significant reduction in autocorrelations, especially in those for the 
topological charge, in terms of wall-clock time used to generate the configurations. 





mia 


m s /mi 


Atught 


T(Ref.[lJ) 


t(MD) 


Acceptance 


(P) 


(VY{m,)) 


V/a 


= 24 3 


x64,L s = 16, p =2.13, 


a~ l = 1.73(3) GeV, 


m res a = 0.003152(43), r/traj = 1 




0.005 


5.3 


1/6 


4460 


8980 


73% 


0.588053(4) 


0.001224(2) 


0.04 




















0.01 


3.3 


1/5 


5020 


8540 


70% 


0.588009(5) 


0.001738(2) 


V/a 


= 32 3 


x 64, L 


= 16, J8 = 2.25, 


a- 1 =2.28(3) GeV, 


m res a = 0.0006664(76), T/traj = 2 




0.004 


6.6 


1/8 




6856 


72% 


0.615587(3) 


0.000673(1) 


0.03 


0.006 


4.6 


1/8 




7650 


76% 


0.615585(3) 


0.000872(1) 




0.008 


3.5 


1/7 




5930 


73% 


0.615571(4) 


0.001066(1) 



TABLE I: Simulation parameters as well as the average acceptance, plaquette ({P}) and value for the light- 
quark chiral condensate ((\j/\j/(mi))) for the ensembles studied in this paper. The fifth column shows the 
number of time units in the ensembles that were included from Ref. [ 1]. The residual masses given explicitly 
and those appearing in the ratio rhi/m s are taken from Table IVIll appearing in Section III below. 



A final optimization was used for the simulations run on the IBM BG/P machines at the Ar- 
gonne Leadership Computing Facility (ALCF). Instead of using double precision throughout, the 
BAGEL- generated assembly routines Ii24ll keep the spin-projected spinors in single precision in 
the conjugate gradient(CG) inverters during the molecular dynamics evolution to decrease the 
amount of communication needed per CG iteration. (Full precision is used in the accept-reject 
step.) While this kind of improvement is expected to make the molecular dynamics integrator un- 
stable for sufficiently large volumes, the effect on the acceptance turned out to be minimal for all 
the ensembles presented in this paper while improving the performance of the CG by up to 20% 
compared to a full double precision CG with the same local volume. 
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FIG. 1 : Evolution of the average plaquette (left panel) and the chiral condensate (right panel) for the /3 = 
2.25, 32 3 x 64, L s = 16 ensembles. The chiral condensate is normalized such that (YY) ~ l/ m i n the heavy 
quark limit. 

B. Ensemble properties 

In Fig. Q] we show the evolution of the plaquette and the chiral condensate for the 32 3 ensembles. 
Both quantities suggest that 500 MD units is enough for the thermalization of each of the 32 3 
ensembles. We have thus begun measurements at 1000 MD units for m/ = 0.006 (except for the 
measurements of the chiral condensate which started after 3304 MD units) and 520 MD units for 
the other 32 3 ensembles. (The starting points for measurements on the three 24 3 x 64 ensembles 
are given in Tab. I of Ref. yj].) 

Figure |2] shows the integrated autocorrelation time for various quantities measured on the 32 3 
ensembles. As can be seen the plaquette, chiral condensate and even the light pion propagator 
for a separation of 20 time units show a short autocorrelation time of 5-10 MD units. However, 
the measured autocorrelation times for the topological charge are much larger, on the order of 
80 MD units. In fact, as is discussed in Section IVlTl the evolutions shown in Fig. [52] suggest 
even longer autocorrelation times implying that the autocorrelation times shown in Fig. [2] may be 
underestimated because of insufficient statistics. 

In Section IVIII this issue of the autocorrelation time for the topological charge is discussed in 
greater detail and the /3 =2.13 and 2.25 evolutions are compared. The 32 3 , /3 = 2.25 ensem- 
bles (with finer lattice spacing) are shown to evolve topology more slowly. This suggests that 
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G-O Plaquette 

B— B Chiral condensate 

A-^. Pseudoscalar(t=20) 

Q—£> Topological charge(m =0.004) 

<j— <| Topological charge(m =0.006) 

<^-^y Topological charge(m =0.008) 




FIG. 2: The integrated autocorrelation time is shown for the average plaquette, chiral condensate (\\fi\f), 
pseudoscalar propagator at time separation 20 from a Gaussian source and point sink, all computed from 
the 32 3 , mi = 0.004 ensemble and the global topological charge for all three 32 ensembles. The chiral 
condensate and plaquette are measured every two MD units and the averages within sequential blocks of 
10 MD units have been analyzed. The topological charge is measured every 4 MD units and the averages 
within sequential blocks of 20 MD units have been analyzed. All other quantities were measured every 20 
MD units and no averaging has been performed. Further discussion of the topological charge is given in 
section IVlTl 



the change from the DBW2 gauge action used in earlier 2-fiavor work I125ll to the Iwasaki gauge 
action used here may have been a wise one. While the DBW2 gauge action gives smaller residual 
DWF chiral symmetry breaking, it does this by suppressing the tunneling which changes topolog- 
ical charge. Thus, the use of the DBW2 gauge action may have resulted in a topological charge 
evolution for our current finest lattice spacings that would have been unacceptably slow. 



C. Fitting procedure 

In the analysis described in this paper it is important to take into account the fact that the various 
quantities computed on a single gauge configuration may be correlated. To do this we apply the 
jackknife technique to simple uncorrelated fits. While there is no proof, or even expectation, that 
this is an optimal procedure, the jackknife will provide a good estimate of the error except in the 
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unlikely event of large deviations of our result from a normal distribution. While we could attempt 
to perform a "text-book" correlated fit (again, using a jackknife procedure), this would not be 
sensible: such fits assume that the data should exactly follow the functional form used in the fit. 
In the case of a fit to chiral perturbation theory or a simpler analytic ansatz for the quark-mass 
dependence of physical quantities we know that this is not the case. While this complaint applies 
to both correlated and uncorrelated fits, for the highly correlated lattice data with which we are 
dealing, small deviations (which in this procedure are assumed to be statistical, but in our case are 
likely to be systematic) are penalized by many orders of magnitude more for the correlated than 
uncorrelated fits. Nevertheless, we have performed correlated fits, where the correlation matrix 
is obtained by taking increasing numbers of the leading eigenvectors. Within our limited ability 
to estimate the correlation matrix, we find no significant difference in the results and errors with 
those obtained using uncorrelated fits. Therefore, in this paper (as was also the case in Ref. yj]) 
we present our main results from the uncorrelated fits, but with a full jackknife procedure for 
estimating the errors. However, it must be borne in mind that for such uncorrelated fits the resulting 
X 2 may not be a reliable indicator of goodness of fit. Therefore, we present a sample set of our fits 
graphically. 

D. Reweighting in the mass of the sea strange-quark 

The sea strange quark mass value used in our ensemble generation, m^ im \ differs from the one 
in nature, which we determine only after performing our final analysis. In this subsection, we 
describe the reweighting method used to correct this strange quark mass from mjf to the tar- 
get mass m^. Various target heavy quark masses are determined in Section IVl through interpola- 
tion/extrapolation to yield meson masses which match either unphysical values present in a dif- 
ferent ensemble or which reproduce those from experiment. Recently, several large-scale QCD 



simulations have been reported using a reweighting technique |26-|28|]. The various uses of this 



method include obtai ning sea quark mass derivatives in Ref. 112 911 . tuning the 



quark masses in Ref. DOh . tuning the strange and charm quark masses in Ref. Bill and going to 



larger L s for the DWF action in Ref. 0211 



ight and strange 



An observable, such as the meson propagator, at the target strange sea quark mass is obtained by 
measuring that observable on the ensemble generated using m^ im \ multiplied by the reweighting 
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factor w: 



(sim) 

m h — / . \ • \ J ) 



Here the reweighting factor w[U^] for a particular ensemble of gauge links is the ratio of the 
square root of the two-flavor Dirac determinant evaluated at the mass m/, divided by that same 
rooted determinant evaluated at mjf , 

det@(m h y/ 2 

w Wu \ — r^—, — r • (6) 

L MJ det^mf 1 ^) 1 / 2 

This factor must be calculated for each configuration on which measurements will be performed 
in the ensemble generated using the sea strange mass m^ im \ 

Among the many possible ways of computing the determinant ratio in Eq. ©, we have chosen to 
use the Hermitian matrix ^(m^m^™- 1 ), whose determinant is wfL^], 



/ (sim) \ 

il(m hl m K h ') 



D(m^ ■ D(m h y ' [D(m h )}-^ D{m* 



1/2 



-1/2 



1/2 



(sim)\ 



1 1/2 



(7) 



The square root of these matrices is implemented using the same rational polynomial approx- 
imation, ^i (x), and multi-shift conjugate gradient algorithm, which are used in the ensemble 

i 

generation. The order of the matrix products in Q. assures that in the limit of — > m^ im \ Q. goes 
to the unit matrix, so that the method described below for evaluating w has vanishing stochastic 
error in this limit. 

To obtain w on each configuration, the determinant of Q. is stochastically evaluated using a com- 
plex random Gaussian vector ^ of dimension L s x 12. Each complex element is drawn from a 
random distribution centered at zero with width in both the real and imaginary directions: 

w _// -r[Q-V(2of)]|vv V 

w-{{e « }}£ = 2 . (8) 

We set ct| = 1/2 and sample using Gaussian vectors per configuration. For one sample, two 
multi-mass inversions, one for m/, and another for mi sl , are performed. 

One needs to be careful in evaluating Eq. © to avoid a large and difficult to estimate statistical 
error. When the eigenvalues of Q., Xq, are far from 1 / (2<j|), the large shift in the width of the 
Gaussian in the integrand will cause poor sampling in this stochastic evaluation of w, as can be 
seen if Eq. [8] is rewritten with Q. diagonal: 

"= EI [ d&fr^-V^^MW) f Y\ jd^{e<^ 2 ° 2 0. (9) 

A^espect(f2) /l^Gspect(fi) 
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The first exponential function in the integrand © will be a rapidly decreasing function of 
when [An — l/(2a|)] is large, with most of the Gaussian samples generated according to the 
second exponential function in Eq. © falling in a region where the first factor is very small. In 
this sense, Eq. ([8]) may provide a statistically noisy estimate of the ratio of the determinants in 
Eq. ©. The fluctuations in this estimate will be rapidly reduced when [Xq — 1/ (2<T?)] — > or, for 
our choice of , when Q. becomes close to the unit matrix, Q. — > 1 . 



To reduce the stochastic noise in our estimate, detQ is divided into N rw factors Il27ll 

Nrw-l Nrw-l 



w = det^= Yl deta i= f] ((g^ tQ, - 1/(2ff « )]§, )) 6 . (10) 

i=0 !=0 

Each of Q.j needs to be close to unit matrix while keeping the determinant of the product the 
same as the original determinant. Each factor dettl, in the product, is evaluated using Eq. © 
with A^s Gaussian vectors. We note that all Gaussian vectors, <fjj, must be statistically independent 
otherwise there will be unwanted correlation among contribution from the Nrw steps. A similar 
decomposition of the reweighting factor is also possible by using the n th root of the operators 02I1 . 
In this work, Q,- is chosen by uniformly dividing the interval [m^, ni^ irc ^\ into smaller pieces: 

Q. i = Q.(ml +l \mf) , (11) 



(sim) 

„(sim) _,_ . m h — m h 

N„ 



K-^h+i Tr h ,(i = 0,l,---,N rw ) . (12) 



In that way, reweighting factors for the intermediate masses are also obtained, which will be 
used in our analysis too. 

For a given difference between the target and the simulation masses, m/, — m^ sim \ Nrw needs to be 
sufficiently large that Q ; is close to the unit matrix, suppressing the statistical noise in estimating 
each of the determinants. We have checked whether N rw is large enough in our calculation of 
the reweighting factor. Figure [3] shows the logarithm of the full reweighting factor, — m(w), as a 
function of the number of divisions in strange quark mass, N rw , on the /3 = 2.13, 24 3 x 64, m; = 
0.005 lattices, the 2,000th trajectory in the left panel and the 4,000th trajectory in the right panel. 
The target and simulation quark masses are = 0.035 and m^ sim) = 0.040. 
For Nrw < 10, the reweighting factor w appears inconsistent with the results obtained for larger N rm 
by a large amount (note that — ln(vv) is plotted) for the left case (2,000th trajectory). We believe 
this is caused by the poor stochastic sampling in our method to compute w when N rw < 10 and 
that for these cases the statistics are insufficient to estimate the error accurately. 
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FIG. 3: Logarithm of the reweighting factor, — ln(w), as a function of the number of divisions in the 
strange quark mass, N rw on the j8 = 2.13 , 24 3 x 64, mi = 0.005 lattices, the 2,000th trajectory on the left 
panel and the 4,000th trajectory on the right panel. The target and simulation quark masses are m/, = 0.035 



and m 



(sim) 



0.040. For N rw = 1,5, 10, 20, 32, 40, the number of Gaussian samples per mass steps is set to 



Ni = 40, 8, 4, 4, 2, 2, respectively. The error bars shown are the standard deviations resulting from N rw x Nt 
samples for det£2,. We interpret the inconsistency between the values for Nrw = 1,5 and 10 and those with 
larger N rw in the left-hand panel as resulting from insufficient statistics leading to under-estimated errors for 
these three cases where the stochastic sampling is very poor. 

ensemble m^ im) m h N rw Nt 



32 3 x 64 0.030 0.025 10 4 
24 3 x 64 0.040 0.030 40 2 



TABLE II: Parameters chosen for the sea strange quark mass reweighting calculation. 

We also check the relative difference between the reweighting factors for Nrw = 20 and N rw = 40 
in Fig.|4]for five lattices. This plot indicates that N rw = 20 is sufficient to estimate the reweighting 
factor and its error for changing from m^ im ^ = 0.040 to m/ 7 = 0.035 on this ensemble. We summa- 
rize the values of N rw and Np used in estimating the reweighting factors for the sea strange quark 
mass in Tab. UH 

Is the N rw dependence, described above, all one needs to check to assure the correctness of the 
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6000 



FIG. 4: The relative differences between the reweighting factors for N rw = 20, N» = 4 and N rw = 40,^ = 2 
on five lattices. The target and simulation quark masses are mh = 0.035 and m£ im) = 0.040. 



reweighting procedure? The answer is clearly no. So far, we have only established that Eq. (1101) 
estimates w to some degree of accuracy, on each configuration for large N rw . One needs further 
checks to see whether or not the reweighted observable in Eq. © has an accurately estimated 
statistical error. A highly inaccurate estimate of the statistical errors could easily result from a 
poor overlap between the reweighted ensemble and the original ensemble generated by the RHMC 
simulation. In addition, because the reweighted observable in Eq. © is given by a ratio of averages 
it is a biased estimator of the observable of interest. In this circumstance, a large statistical error, 
even if well determined, may lead to a systematic error of order l/Af con f enhanced by this large 
statistical error. 

We have attempted the following checks: In Fig. [51 w is plotted as a function of trajectory. If the 
fluctuation among different configurations is large, Eq. © might be dominated by a small number 
of measurements made on those configurations with large w, and the measurement efficiency for 
the reweighted observable would be very poor. Using the reweighting factor, w;, obtained on the 
i ,h configuration, the reweighted observable & can be written from Eq. © as, 

Nconf 



I 0i 



Wj 



Wi 



i=l 
Wi 



(13) 
(14) 



Because the process of reweighing selectively samples the original distribution, even with pre- 
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ensemble 


max(w^ 


min(w^ 




^onf 


24 3 


x64,m, = 0.005 


10.0 


0.078 


90.3 20.3 


203 


24 3 


x 64, mi = 0.010 


5.50 


0.049 


97.0 32.4 


178 


32 3 


x 64, m, =0.004 


4.77 


0.17 


228 63.9 


305 


32 3 


x 64, mi =0.006 


3.45 


0.23 


234 90.4 


312 


32 3 


x 64, m, =0.008 


5.36 


0.16 


183 47.0 


252 



TABLE III: The maximum and minimum reweighting factors, the effective number of samples, N e ff, ac- 
cording to the formula derived in this paper, (Eq. ([T5l)). the corresponding number, N* s given by the 



formula of Ref. H33TI (defined in Eq. (fT6l >) and the actual number of configurations -/V CO nf in each en- 



semble. The target sea strange quark mass and that of the simulation are mn = 0.0345, mf m) = 0.040 
(m h = 0.0275, mf m) = 0.030) for 24 3 x 64 (32 3 x 64). 

cisely determined reweighting factors we should expect the effective number of samples to be 
reduced and the statistical errors to increase. In Appendix O this effect is analyzed in the case that 
correlations between the data and the reweighting factors can be neglected when estimating these 
statistical errors, including the effects of autocorrelations. For the case of no autocorrelations, we 
obtain the following expression for the effective number of configurations after reweighting: 

Neff= V „ - / ■ (15) 

The quantity A^ff goes to Af con f if there is no fluctuation in the w; while it goes to 1 if the largest 
Wi completely dominates the reweighted ensemble. We summarize the statistical features of the 
reweighting factors for each ensemble in Tab.lTTTl For completeness we also compare the definition 



of A/eff given in Eq. (1151) with the more pessimistic estimate used in Ref. [13311 : 



Kff= f"\ ■ (16) 
max j[wj) 

As can be seen from Tab. [Till our choice gives a somewhat more optimistic view of the effects of 
reweighting on the effective size of our ensembles. 

As the numbers in Tab. ITTT1 indicate, for our ensemble and reweighting settings, the ensembles are 
not overwhelmed by a small number of configurations. 

The efficiency of the reweighting procedure is also observable dependent. It is influenced by the 
fluctuations of the reweighted observable within the ensemble and the strength of the correlation 
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w at m. =0.0345 from 0.040 

norm h 

m,=0.005, 0.010 




w at n\ =0.0275 from 0.030 

norm h 

m.=0.004, 0.006, 0.008 



1 
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FIG. 5: The normalized reweighting factor w; as a function of trajectory number i for the 24 3 x 64, m/ = 
0.005, 0.010 ensembles (left-hand plot) and the 32 3 x 64, m, = 0.004, 0.006, 0.008 ensembles (right-hand 
plot). The sea quark masses ni/ are plotted in ascending order from top to bottom. The target sea strange 
quark mass and that of simulation are m h = 0.0345, mf m) = 0.040 (m h = 0.0275, mf m) = 0.030) for the 
left-hand (right-hand) plot. 



between the reweighted observable and the reweighting factor. Sanity checks of the statistical 
properties of the most important observables, m % and f % , have been performed and are summarized 
in Fig. [6l The observables reweighted to w/, = 0.0250 from m^ im ^ = 0.030 are calculated using 
the first half and the second half of the ensemble (circle symbols), which are compared to that 
of the full statistics (square symbols). The number of the Gaussian vectors, Np, is also varied 
from Np = 1 (blue symbols) to Np =4 (red symbols) in the same plot. In the case of m K , all the 
statistical samples are within 1 x o, while for f K the deviations are less than ~2xrj. 
To probe the /??/, dependence of the observables, we show in Fig. [7] the correctly reweighted m K 
and f n as a function of m/, along with the results obtained from randomly permuting the {w,} in 
Eq. (fT3l) . The random permutation is done for each reweighted mass m/ 7 to show the difference 
from the correctly reweighted observables. While the randomly reweighted observables are almost 
flat in m/,, the correctly reweighted observables have a positive slope in m/,. Finally in Fig. [8] we 
plot the reweighted observables f n and as a function of the target reweighted mass for three 
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Unitary m versus hit index reweighted to 0.02500 



Unitary f versus hit index reweighted to 0.02500 
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FIG. 6: Reweighted values for m n (left) and /jt (right) for various numbers of reweighting hits, = 
1 (blue), Np =2 (green), Np = 4 (red) ) on each ensemble. The squares are for the full data set (300 
configurations) and the circles are for the first and second half of the data (150 configurations.) The data is 
from the 32 3 x 64 x 16, (m/,m/j) = (0.004,0.03) ensemble with a light valence quark of mass 0.004. The 
black symbols are the unreweighted observables. 
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FIG. 7: The left figure gives m % with correct reweighting factors (blue squares) and with randomly permuted 
reweighting factors (green diamonds). The right figure is the same but for f n . 

example parameter points. Note that in both Figs. [7] and [8] we observe an increase in statistical 
errors which appears roughly consistent with what should be expected from the decrease in vWeff- 
We should emphasize that further careful studies may be needed to establish a more accurate 
estimate of possible errors in the reweighting procedure. 
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Unitary f versus m h from reweighting 



Unitary f versus m h from reweighting 
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FIG. 8: Reweighted results for f n (left) and f% (right) as functions of ni/, at three parameter sets (/$,m/): 
green diamonds: (2.25, 0.008), red circles: (2.13, 0.005), blue squares: (2.25, 0.004). 

Volume (m/,m/,) Total MD time Measurement range Measurement total 



24 3 


(0.005, 0.04) 


0-8980 


900-8980 every 40 


203 


24 3 


(0.01,0.04) 


1455-8540 


1460-8540 every 40 


178 


32 3 


(0.004, 0.03) 


0-6756 


520-6600 every 20 


305 


32 3 


(0.006, 0.03) 


0-7220 


1000-7220 every 20 


312 


32 3 


(0.008, 0.03) 


0-5930 


520-5540 every 20 


252 



TABLE IV: Summary of the five ensembles used in this work. 



III. UPDATED RESULTS FROM THE 24 3 ENSEMBLES 



In this section we update the results presented on the 24 3 ensembles in yj] to the extended data set 
described in Sec. [Til and in TableUin particular. For this extended data set we make measurements 
of pseudoscalar quantities on a total of 203 configurations for the m/ = 0.005 ensemble and 178 
configurations for the m/ = 0.01 ensemble. These configurations were separated by 40 trajecto- 
ries as documented in the first two rows of TablelTVl In our previous work we used 92 of these 
measurements on each ensemble III, |4J]. Before performing the analyses we binned the data into 
blocks of either 80 or 400 trajectories and the measurements from each bin were then treated as 
being statistically independent. No statistically significant increase in the error was observed with 
the analysis using bins of 400 trajectories compared to that with bins of 80 trajectories. 



19 



In the following sections the results from the 24 3 lattices, combined with those obtained on the 
32 3 ensembles, will be input into global chiral and continuum fits in order to determine physical 
quantities; here we simply tabulate the fitted pseudoscalar masses and decay constants as obtained 
directly from the correlation functions at our simulated quark masses. In addition, since we use 
the mass of the Q. baryon in the definition of the scaling trajectory, we also present the results for 
rrikhh here together with those for the Sommer scale tq and also the scale r\. Finally, in Sec MI Al 
we give the results for the masses of the nucleons and A baryons from the 24 3 ensembles, although 
the chiral and scaling behaviour of these masses will not be studied in this paper. We present these 
baryon masses partly for completeness and partly to share our experience in the use of different 
sources. 

On the 24 3 lattices discussed in this section, the measurements are presented for the two values 
of the sea light-quark mass, m/ = 0.005 and 0.01, and for the full range of valence quark masses 
m = 0.001, 0.005, 0.01, 0.02, 0.03 and 0.04. The ensembles with m, = 0.02 and 0.03, presented 
in yj], are not included in this paper because such values of m/ were found to be too large for 
SU(2) chiral perturbation theory to describe our data. The value of the sea strange-quark mass in 
these simulations is m/ 7 = 0.04. After completing the global chiral and continuum fits described 
in Section|V] below, we find that the physical value of the bare strange-quark mass, obtained using 
the chiral perturbation theory ansatz, is m s = 0.0348(11). In this section we anticipate this result 
and use reweighting to obtain results also at this value of the strange-quark mass. 
For the 24 3 ensembles, we placed Coulomb gauge-fixed wall sources at t = 5 and at t = 57. From 
each source, we calculated two quark propagators, one with periodic and the other with anti- 
periodic boundary conditions. From the periodic propagators for the two sources, denoted by Dp^ 
and Dp^-j, and the anti-periodic propagators, written as D^ l 5 and D A l 51 , we form the combinations 

D KAj = l{ D V+ D £) and D ^57 = 5(^57+^,57) • ( 17 ) 

The use of periodic plus anti-periodic boundary conditions in the time direction doubles the length 
of the lattice in time, which markedly reduces the contamination from around-the- world propaga- 
tion in the time direction. For two point functions, such as the propagator of a pseudoscalar meson 
given by 

(^(0^(0))=£Tr|[D p | A5 (?,x)] t D p ; A5 (r,x)| , (18) 
on a lattice of time extent N t the time dependence of the contribution of the ground state is given 
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by 

(n(t)n(O)) = A [exp(-m K (t - 5)) + exp(-m 7r (2JV* - (t - 5))] . (19) 

Here A is a f-independent constant. For our 24 3 ensembles, we find that around-the-world propa- 
gation is not visible in two-point functions. This is not the case however, for three-point functions, 
as we now explain (although we do not analyze three-point functions in this paper, they are being 
evaluated in the computation of Bk, for example [34]). 

For three-point functions of the form (P(x)0(y)P(z)), where P(x) and P(z) are pseudoscalar in- 
terpolating fields and 0(y) is an operator whose matrix element we wish to measure, we use the 
wall source at t = 5 as the source for P(z) and the wall source at t = 57 as the source for P(x). We 
only consider yo in the range 5 < yo < 57, so we do not perform any measurements in the doubled 
lattice. The doubling of the lattice is important to reject around-the-world propagation in time for 
such measurements. For kaons, we found that a time separation of 52 between the sources gave us 
a broad plateau, with sufficiently small errors. This measurement strategy was chosen to optimise 
the measurement of the kaon bag parameter J4 , 34 ] . 



Before presenting our results for masses, decay constants and r$ and r\, we discuss the values of 
the residual mass and the renormalization constant of the local axial current. The residual mass 
m res( m /) at eacri partially quenched valence mass used in this work is measured using the ratio 117 911 

mJmf) = (W' (20) 
where 7| ? is the usual DWF mid-point pseudoscalar density composed of fields of each chirality 
straddling the mid-point in the fifth dimension, and /| is the physical pseudoscalar density at the 
surfaces of the fifth dimension composed of surface fields in the fifth dimension. The results are 
given in Table |VJ For completeness we also present the corresponding residual masses obtained 
after reweighting to the physical strange mass in Table [VTl The residual mass in the two-flavor 
chiral limit m res = m' KS (m x = mi — 0) is given in Table lVUl and in the left-hand plot of Figure [9l 
We define Za to be the renormalization constant of the local axial current, A^, composed of the 
physical surface fields. Here we have determined Za through two methods. In the first, is 
determined for each valence mass using the improved ratio Il35n of the matrix element {s^4(t)P(0)) 
to (A<x(t)P(Q)), where £^ is the conserved DWF axial current and the results are presented in 
Table I Villi This method assumes = 1, and we find Z& = 0.71651(46) in the two-flavor chiral 
limit with the simulated sea strange mass, and Za = 0.71689(51) when reweighted to the nearby 
physical strange mass. This determination of Za is illustrated in the plots of FigureQIjl As pointed 
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m x 


mi 




0.005 0.01 


0.001 


0.003194(16) 0.003286(28) 


0.005 


0.003154(15) 0.003259(26) 


0.01 


0.003079(14) 0.003187(24) 


0.02 


0.002939(12) 0.003042(21) 


0.03 


0.002822(12) 0.002919(19) 


0.04 


0.002725(11) 0.002818(17) 



TABLE V: m' res (m x ) measured on the 24 3 ensembles at the simulated strange quark mass m/, = 0.04. 



m x 


m t 




0.005 0.001 


0.001 


0.003146(27) 0.003224(33) 


0.005 


0.003099(27) 0.003191(32) 


0.01 


0.003025(26) 0.003120(31) 


0.02 


0.002889(24) 0.002981(26) 


0.03 


0.002774(23) 0.002863(23) 


0.04 


0.002680(21) 0.002765(21) 



TABLE VI: m' res (m x ) on the 24 3 ensembles at the physical strange quark mass. 

out in yj], we expect Z^ = 1 + 0(am ies ), and in yj] we added a ~ 1% error to account for the 
size of this correction. As part of our current work, we have investigated the consequences of 
this correction, which is discussed in detail in appendixE From this analysis, we find Za = 
0.7041 (34), a 1.8% difference from the result with our previous method. Although, as we will see, 
this error is smaller than our current combined errors on the decay constants and other physical 
quantities, we choose to use this value of Za = 0.7019(26), coming from Zy/Zy as defined in 
Equation ( IB 191) , as the normalization factor for the local axial current when quoting all our central 
values below. Here V and Y arc the local and conserved vector currents. 

We now turn to the measurements of the meson masses and decay constants. In order to illustrate 
the quality of the fits, we start by presenting some sample plots for the unitary pion and kaon on 
the in/ = 0.005, = 0.04 ensemble. The pion effective masses obtained using different sources 
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m h 



ml 



in 



phys 



in 



24" 



in 



32 3 



0.003152(43) 0.0006664(76) 
0.003076(58) 0.0006643(82) 

TABLE VII: m res in the two-flavor chiral limit on the 24 3 and 32 3 ensembles at the simulated and physical 
strange sea-quark masses. 
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FIG. 9: Chiral extrapolation of the unitary values of m' Tes for the 24 3 (left) and 32 3 (right) ensembles. 
While the fit is only marginally acceptable for the 32 3 lattices, an additional uncertainty of 0(5 x 10~ 6 ) is 
negligible. 

and sinks are shown in Figure[HJ The mass and decay constant is obtained from a simultaneous fit 
with a single, constrained mass to five correlation functions. These are the (P\P), (A\A) and (A\P) 
correlation functions (denoted in the figure by PP, AA and AP respectively) with gauge-fixed wall 
sources and local (LW) or wall (WW) sinks (we do not use the AA-WW combination because 
it is noisier). The long time extent N t = 64 on our lattices together with the noise properties of 
pseudoscalar states allow for long plateaux and the results are insensitive to the choice of t m m , the 
starting point of the fits. Figure[[2]displays the effective masses for the unitary kaon, together with 
the results obtained from a simultaneous constrained fit. We give an example of the mh dependence 
of the unitary pion and kaon masses in figure [T3J This dependence is obtained by reweighting. 
We normalize the states so that, for periodic boundary conditions, the time dependence of the 
correlators for large times is given by 



(OjO* 1 j7r><7rjO* 2 |0> 



IrrixyV 



e -m xy t _|_ e -m xy {2N,-t) 



(21) 



where the superscripts specify the type of smearing and the subscripts denote the interpolating 
operators. The sign in the square brackets in Eq. (l2~TT) is + for PP and AA correlators and — for AP 
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Z A (chiral) 


Z A (mi =0.005) 


Z A ( mi = 0.01) 


mf m = 0.04 


0.71651(46) 


0.71732(14) 


0.71783(15) 


phys 


0.71689(51) 


0.71746(17) 


0.71781(17) 



TABLE VIII: Z A on the 24 3 ensembles at the simulated and physical strange sea-quark masses. 
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FIG. 10: Measurement of Z A for m/ = 0.005 on the mi = 0.005, m/, = 0.04 ensemble (left panel) and the 
unitary chiral extrapolation of Z A for the 24 3 ensembles (right panel). The results do not change significantly 
under reweighting to the physical strange mass. 



ones. We therefore define the amplitude of the correlator to be 



_ (0\O?\k){k\(%\0) 



IrrixyV 



(22) 



For each correlator included in the simultaneous fit 



J/S^^S^ and 



A/WW 

^ Y AP ) 



we determine the amplitude and obtain the decay constant f xy using 



m xy jV™ 



(23) 



TablellXlcontains the measured pseudoscalar masses and decay constants at the simulated strange- 
quark mass m/, = 0.04. After reweighting to the estimated physical strange-quark mass m s = 
0.0348(1 1) the masses and decay constants of the pions are presented in TablelXland those for the 
kaons in TablelXll 

The Q. baryon, being one of the quantities included in the definition of our scaling trajectory (see 
Section|V]), plays an important role in our analysis. We have performed measurements on the same 
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FIG. 1 1 : Effective pion masses from the PP LW correlator (top left), PP WW correlator (top right), AP LW 
correlator (center left), AP WW (center right) and AA LW correlator (bottom). Note the different vertical 
scale for the WW correlators. The horizontal bands represent the result for the mass from a simultaneous 
fit. 



configurations using a gauge-fixed box source of size 16 lattice units that gives a good plateau 
for the Estate for valence quark masses m x = 0.04 and m x = 0.03 to enable interpolation to the 
physical strange-quark mass. We display the fit to the m x = 0.04 £1 baryon mass on the mi = 0.005, 
m/j = 0.04 ensemble in figureO along with the dependence of this mass on the dynamical strange 
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FIG. 12: Effective kaon masses from the PP LW correlator (top left), PP WW correlator (top right), AP LW 
correlator (center left), AP WW (center right) and AA LW correlator (bottom). Note the different vertical 
scale for the WW correlators. The horizontal bands represent the result for the mass from a simultaneous 
fit. 



mass using reweighting. 

The results for the Q. mass, mhhh, obtained directly at the simulated strange-quark mass (m/, = 0.04) 
with valence strange-quark masses m y = 0.04 and 0.03 are presented in Table lXIIl In this table we 
also present the results for ln^h obtained after reweighting to the physical strange-quark mass. In 
Table lXIlH we display the values of the Sommer scale rrj, r\ and their ratio at both the simulated 
and physical strange-quark masses. These quantities were determined using Wilson loops formed 
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FIG. 13: We illustrate the dependence of the unitary pion (left panel) and kaon (right panel) masses 
on the mi = 0.005, 24 3 ensemble. The values are obtained by reweighting around the simulated value 
(m h = 0.04). 
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FIG. 14: Fit to the Q. baryon mass with valence strange mass m x = 0.04 on the m, = 0.005, m h = 0.04, 24 3 
ensemble showing the quality of the fit with our box source (left panel). We also show the weak dependence 
of the Q. baryon mass with fixed valence mass m x = 0.04 on our simulated mu inferred by the reweighting 
procedure on the mi = 0.005, 24 3 ensemble (right panel). 

from products of temporal gauge links with Coulomb gauge-fixed closures in spatial directions, 
with an exponential fit to the time-dependence of the Wilson loop W(r,t) from t = 3 to t = 7 for 
each value of the separation r. The resulting potential V(r) was then fit over the range r = 2.45 — 8 
to the Cornell form Bol 

(X 

V(r)=V -- + ar, (24) 
r 

where Vq, oc and o are constants. These fits are illustrated in Figure [151 which shows the fit to the 
time dependence of the Wilson loop W(r = 2.45,?) at the physical strange-quark mass, and also 
the subsequent fit over the potential. The strange-quark mass dependence of the scales ro and r\ is 
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m r m v 



0.04 0.04 

0.03 0.04 

0.02 0.04 

0.01 0.04 

0.005 0.04 

0.001 0.04 

0.03 0.03 

0.02 0.03 

0.01 0.03 

0.005 0.03 

0.001 0.03 

0.02 0.02 

0.01 0.02 

0.005 0.02 

0.001 0.02 

0.01 0.01 

0.005 0.01 

0.001 0.01 

0.005 0.005 

0.001 0.005 

0.001 0.001 



(0.005) m„(0.01 



0.4317(4 
0.4051(4 
0.3772(5 
0.3478(5 
0.3325(6 
0.3199(7 
0.3771(4 
0.3472(5 
0.3152(5 
0.2983(5 
0.2843(6 
0.3149(5 
0.2794(5 
0.2603(5 
0.2440(6 
0.2389(5 
0.2161(5 
0.1960(6 
0.1904(6 
0.1669(6 
0.1391(6 



0.4344(4 
0.4080(4 
0.3802(4 
0.3509(5 
0.3358(5 
0.3233(7 
0.3800(4 
0.3502(4 
0.3184(4 
0.3016(5 
0.2877(6 
0.3179(4 
0.2826(5 
0.2636(5 
0.2475(6 
0.2422(5 
0.2195(5 
0.1997(6 
0.1940(6 
0.1709(6 
0.1434(7 



.£,,(0.005) .U0.01) 



0.1063(6) 0.1087(6) 

0.1034(6) 0.1059(6) 

0.1002(5) 0.1028(5) 

0.0967(5) 0.0996(6) 

0.0949(5) 0.0982(6) 

0.0937(6) 0.0975(7) 

0.1006(5) 0.1031(5) 

0.0974(5) 0.1001(5) 

0.0939(5) 0.0969(5) 

0.0920(5) 0.0954(6) 

0.0908(6) 0.0946(6) 

0.0943(5) 0.0971(5) 

0.0908(5) 0.0938(5) 

0.0889(5) 0.0923(5) 

0.0876(5) 0.0915(6) 

0.0872(5) 0.0905(5) 

0.0853(5) 0.0889(5) 

0.0840(5) 0.0879(5) 

0.0834(5) 0.0871(5) 

0.0819(5) 0.0858(5) 

0.0802(5) 0.0840(5) 



TABLE IX: Pseudoscalar masses m^mi) and decay constants fxy(mi) on the 24 3 ensembles at the simulated 
strange-quark mass (m/, = 0.04). 



small and cannot be resolved within our statistics. 



A. Nucleon and A Masses 



A detailed study of the baryon mass spectrum, including the continuum and chiral extrapolations, 
is postponed to a separate paper. The one exception is the Q. baryon, whose mass is used in the 
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m x m y 


m^O.005) m^(O.Ol) 


^(0.005) /^(O.Ol) 


0.01 0.01 


0.2378(8) 0.2420(7) 


0.0867(5) 0.0900(6) 


0.005 0.01 


0.2149(9) 0.2192(7) 


0.0848(6) 0.0882(6) 


0.001 0.01 


0.1948(10) 0.1994(8) 


0.0833(6) 0.0871(6) 


0.005 0.005 


0.1891(10) 0.1936(8) 


0.0828(5) 0.0863(6) 


0.001 0.005 


0.1656(11) 0.1704(8) 


0.0813(6) 0.0850(6) 


0.001 0.001 


0.1377(12) 0.1427(9) 


0.0796(6) 0.0832(7) 



TABLE X: Pion masses m^mi) and decay constants fxy(mi) on the 24 3 ensembles at the physical strange- 
quark mass m s = 0.0348(11). 



m x 


m xh (0.005) m^(O.Ol) 


fxh (0.005) 


MOM) 


0.01 


0.330(4) 


0.334(4) 


0.0947(7) 


0.0978(8) 


0.005 


0.314(4) 


0.318(4) 


0.0928(7) 


0.0963(9) 


0.001 


0.301(4) 


0.305(4) 


0.0915(8) 0.0955(10) 



TABLE XI: Kaon masses m x i t (mi ) and decay constants f x i, (mi ) on the 24 3 ensembles at the physical strange- 
quark mass m s = 0.0348(11). 

definition of the scaling trajectory and which is therefore studied in detail together with the prop- 
erties of pseudoscalar mesons. In this subsection we briefly discuss our experiences in extracting 
the masses of the nucleons and A-baryons using different sources and present the results for these 
masses on each ensemble, starting here with those from the 24 3 ensembles. The baryon spec- 
trum from the 32 3 ensembles will be discussed in Sec. lIV Al We start however, with some general 
comments about our procedures which are relevant to both sets of ensembles. 
We use the standard operator, N = e a b c (u T a Cy^db)u c , to create and annihilate nucleon states and 
A = £ a bc( u aCYv- u b) u c for the flavor decuplet A states. On an anti -periodic lattice of size N t in the 
time direction, the zero-momentum two-point correlation function, C(t), calculated with one of 
these baryonic operators at its source and sink, takes the following asymptotic form for sufficiently 
large time, t, 

C(t) = Z[(l + YA)e- Mt - (1 - Y4)e- M{N '-% (25) 

corresponding to particle and anti-particle propagation, respectively. Conventionally one chooses 
an appropriate range in time where the excited-state contributions can be neglected so that this 
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m y 


m h 


m a (0.005) m n (0.01) 


0.04 


0.04 


1.013(3) 1.028(4) 


0.03 


0.04 


0.963(4) 0.978(4) 


0.0348 0.0348 


0.988(9) 1.001(7) 



TABLE XII: Omega baryon masses on the 24 3 ensembles at the simulated strange quark mass m/, = 0.04 
(first two rows) and at the physical strange-quark mass (third row). 



Quantity 



m h = 0.04 
2(0.005) 2(0.01) 



m h = 0.0348 
2(0.005) 2(0.01) 



r 4.16(2) 4.10(2) 4.15(2) 4.12(3) 
2.82(3) 2.70(2) 2.83(3) 2.72(3) 

ri/r 

TABLE XIII: The quantities ro, r\ and n/ro at the simulated (m/, = 0.04) and physical (m/, = 0.0348) 
strange quark masses on the 24 3 ensembles. Q(mi) denotes the quantity measured with light-quark mass 
mi. 

form is valid, and extracts the ground-state mass, M, by fitting the numerical data to the function 
in Eq. (|25l) . This is indeed what we do to extract baryon masses from the 24 3 ensembles. Alter- 
natively we can try to fit the correlation function to a sum of two exponentials, representing the 
ground- and excited-state contributions. As will be reported in Sec. lIV Al this is the method we 
use for the 32 3 ensembles. 

The determination of baryon masses can be made more effective by an appropriate choice of 
smearing at the source and/or sink. We use several different choices of the smearing of these 



operators, wall, box, and gauge-invariant Gaussian H37L I38H . m an attempt to obtain a better overlap 
with the ground state; our choices are summarized in Table IXIVI The wall source, used for the 
32 3 ensembles, is Coulomb-gauge fixed. A box source of size 16, also Coulomb-gauge-fixed, is 
used for the 24 3 ensembles. The Gaussian-source radius is set to 7 lattice units and 100 smearing 
steps are used for the 24 3 ensembles, while the radius is 6 in the 32 3 ensembles: these choices are 
optimized for our nucleon-structure calculations B7H9J]. 

As can also be seen from the table, several steps are taken to reduce the statistical error. For each 
configuration, as many as four different time slices are used for the sources, usually separated 
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o.5 - e 



FIG. 15: The effective potential of the Wilson loops with a spatial extent of r = 2.45 on the 24 3 , mi = 0.005 
ensemble at the physical strange-quark mass, overlaid by the fit to the range t = 3 — 7 (left panel). The right 
panel shows the static inter-quark potential V (r) on this ensemble, again at the physical strange-quark mass, 
as a function of the spatial extent of the Wilson loops, overlaid by the fit to the Cornell form over the range 
r = 2.45-8. 

by 16 lattice units, but occasionally fewer. Measurements are made as frequently as every tenth 
trajectory and are averaged into bins of 40 hybrid Monte Carlo time units. 

We now turn to the results obtained specifically on the 24 3 ensembles. The unitary nucleon and 
A effective masses are plotted in Figs. [16] and [FT] for each choice of quark mass. For the nucleon, 
both Gaussian and box sources are shown. Plateaus for the effective masses obtained with the 
box source appear quickly, suggesting a strong overlap with the ground state. The corresponding 
plateaus obtained with the Gaussian source appear more slowly, from above. Both sets of results 
agree reasonably well for sufficiently large t. For the A the correlators were only computed using 
the box source and the plateaus for the effective masses again appear quickly. The results for 
the masses, obtained using fully correlated fits, are summarized in Table IXVI Note such fully 
correlated fits work well for extracting baryon masses as the procedure involves much shorter 
ranges in time than for the meson observables discussed in the rest of this paper. As expected from 
the effective mass plots, nucleon masses obtained using different sources agree fairly well when 
the fits are performed over appropriate ranges. All values of % 2 ld.o.f. are close to 1 or smaller, 
except for the box-source nucleon fit at mj = 0.02 which is about 2.5. 



Some of these results have been reported earlier at Lattice 2008 [390, and also partially in related 
papers on nucleon structure [8,9]]. A preliminary report on a bootstrap correlated analysis with 



frozen correlation matrix was presented at Lattice 2009 ll40ll and the results agree with the updated 
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size mi source type correlators source time slices configurations 



24 3 


0.005 Gaussian 


N 


0,8,16,19,32,40,48,51 647 




0.005 Box 


A, £2 


0,32 90 




0.01 


Gaussian 


N 


0,8,16,19,32,40,48,51 357 




0.01 


Box 


a, a 


0,32 90 




0.02 


Gaussian 


N 


0,8,16,19,32,40,48,51 99 




0.02 


Box 


A, £2 


0,32 43 




0.03 


Gaussian 


N 


0,8,16,19,32,40,48,51 106 




0.03 


Box 


a, a 


0,32 44 


32 3 


0.004 Gaussian 


N, A 


10, 26, 42, 58 264 




0.004 Wall 


N, A 


0, 16, 32, 48 305 




0.006 Wall 


N, A 


0,16,32, 48 224 




0.008 Gaussian 


N, A 


10, 26, 42, 58 169 




0.008 Wall 


N, A 


0, 16, 32, 48 254 



TABLE XIV: Summary of the configurations used in the calculation of the baryon spectrum. 

mi N (Gaussian) N (Box) A (Box) 

0.005 0.671(4) {6-12} 0.669(7) {4-12} 0.865(11) {4-12} 
0.01 0.699(5) {9-15} 0.706(6) {4-12} 0.891(8) {4-12} 
0.02 0.800(8) {8-15} 0.803(7) {4-12} 0.963(8) {4-12} 
0.03 0.896(7) {8-15} 0.894(8) {5-12} 1.029(12) {5-12} 

TABLE XV: Baryon mass in lattice units from the /3 = 2.13, 24 3 ensembles. {} denotes fit range, 
ones given here. 

IV. RESULTS FROM THE 32 3 ENSEMBLES 

The results for masses, decay constants, ro and r\ obtained directly on the 32 3 lattice are pre- 
sented in the same format as those from the 24 3 ensembles in Section!]]] and the available 
measurements are also summarised in table [IV] The results are presented for three values of 
the sea light-quark mass m/ = 0.004, 0.006 and 0.008 which correspond to unitary pion masses 
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FIG. 16: Nucleon effective mass plots from the 24 3 ensembles. Results obtained using the Gaussian source 
are marked by red squares and those from the box source by blue circles. The four plots correspond to 
unitary light-quark masses 0.005 (top-left), 0.01 (top-right), 0.02 (bottom-left) and 0.03 (bottom-right). 

in the range 290 MeV - 400 MeV which we had found to be consistent with SU(2) chiral per- 
turbation theory on the 24 3 lattice yj]. The valence-quark masses used in the analysis are 
m Xjy = 0.002, 0.004, 0.006, 0.008, 0.025 and 0.03. For pseudoscalar quantities we use 305, 312 
and 252 measurements separated by 20 trajectories on the 0.004, 0.006 and 0.008 ensembles re- 
spectively (see TablelTvT). For the 32 3 lattices, we have used a single-source technique for our 
measurements of pseudoscalar quantities, which differs from the two-source method for the 24 3 
ensembles. Recall that for the 24 3 ensembles, as discussed in Section!]]]], we placed Coulomb 
gauge-fixed wall sources at t = 5 and at t = 57. For the 32 3 ensembles we have used a sin- 
gle source and calculated both periodic and anti-periodic propagators from this one source. The 
source is placed at t = on the first configuration used for measurements, and the position of 
the source is then increased by 16 for every subsequent measurement so that f src = \6n mod 64 
where n is the measurement index, which starts from zero. Moving the source in this way helps 
to decorrelate measurements. We always place the anti-periodic boundary condition on the links 
in the time direction going from the hyperplane with t = t src — 1 to t = t SIC . Clearly the number of 
propagators to calculate for the single source method is half that for the two-source method. 
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FIG. 17: Effective mass plots for the A baryon from the 24 3 ensembles. The results were obtained using the 
box source. The four plots correspond to unitary light-quark masses 0.005 (top-left), 0.01 (top-right), 0.02 
(bottom-left) and 0.03 (bottom-right). 

For meson two-point functions, as given in Eq. (fT8l) , the single-source method is identical to the 
two-source method, except for having half the number of measurements per configuration. For the 
light-quark masses on our 32 3 ensembles we do see around-the-world effects at the fraction of a 
percent level, so fits of the form in Eq. ( fT9l must be used. We also perform measurements using 
three-point functions of the type (P(x)0(y)P(z)), where P(x) and P(z) are pseudoscalar interpo- 
lating fields and 0(y) is an operator whose matrix element we wish to measure. Here P(x) is made 
out of propagators of the form D p ^ A Q = 1/2 ^J+D^qJ in the notation of Eq. (fT71) and P(z) 
is composed of Dp\ A = 1/2 (d^q —D a 1 qJ propagators. This means that the time separation be- 
tween P(x) and P(z) is N t , the time extent of our lattice. We performed tests on our 24 3 ensembles, 
comparing the single-source and two-source methods and found that, for the same number of in- 
versions, the single-source methods gave at least as small an error as the two-source methods. The 
single-source method allows us to measure on more configurations for the same computer time 
and so we chose this method. Although we do not discuss three-point measurements in this paper, 
sharing propagators between them and the two-point measurements discussed here has helped to 
define our measurement strategy. 
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m x 


mi 




0.004 0.006 0.008 


0.002 


0.0006761(35) 0.0006688(34) 0.0006822(37) 


0.004 


0.0006697(34) 0.0006651(31) 0.0006791(36) 


0.006 


0.0006622(33) 0.0006589(30) 0.0006736(35) 


0.008 


0.0006550(32) 0.0006524(29) 0.0006676(34) 


0.025 


0.0006090(24) 0.0006089(21) 0.0006218(25) 


0.03 


0.0005993(23) 0.0005997(20) 0.0006115(24) 



TABLE XVI: m' res on the 32 3 ensemble set at the simulated strange quark mass my, = 0.03. 



m x 


mi 




0.004 0.006 0.008 


0.002 


0.0006718(39) 0.0006671(36) 0.0006781(44) 


0.004 


0.0006658(39) 0.0006633(33) 0.0006751(42) 


0.006 


0.0006586(37) 0.0006569(31) 0.0006696(40) 


0.008 


0.0006515(36) 0.0006503(30) 0.0006636(39) 


0.025 


0.0006063(26) 0.0006058(24) 0.0006180(31) 


0.03 


0.0005967(24) 0.0005966(22) 0.0006080(29) 



TABLE XVII: m' res on the 32 3 ensemble set at the physical strange quark mass. 

The measured values of the residual mass m' ies at each pair of valence and sea light-quark masses 
(m x ,m[) used in this work are given in table IXVlt in this table the strange-quark mass is the 
one used in the simulation m/, = 0.03. Table IXVTIl contains the corresponding results obtained 
after reweighting to the physical strange mass (m s = 0.0273(7)) determined later in the analysis 
and presented in Section[V] The residual mass in the unitary two-flavor chiral limit is given in 
table EH and figure gj 

The results for Za for the 32 3 ensembles obtained from the ratios of matrix elements of g/4 and A4 
are given in table IXVlTIl We obtain Za = 0.74475(12) in the chiral limit with the simulated sea 
strange mass and Za = 0.74468(13) when reweighted to the nearby physical strange mass. This is 
illustrated in figure [[8] As explained in Section|ni]and appendix|B] however, in this paper we use 
Zy jZy = 0.7396(17) as the normalization factor for the local axial current when calculating the 
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in), 



0.03 



in 



phys 



Z A (chiral) Z A {m x = 0.004) Z A {m x = 0.006) Z A (nn = 0.008) 



0.74475(12) 0.745053(54) 0.745222(45) 0.745328(48) 
0.74469(13) 0.745059(52) 0.745239(47) 0.745384(56) 



TABLE XVIII: Z A on the 32 3 ensembles at the simulated and physical strange sea-quark masses. 
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FIG. 18: Measurement of Z A for my = 0.004 on the m/ = 0.004, m/, = 0.03 ensemble (left panel) and 
the unitary chiral extrapolation of Z A for the 32 3 ensemble set (right panel). The results do not change 
significantly under reweighting to the physical strange mass. 

central values of physical quantities. 

In order to illustrate the quality of the fits, we present sample effective mass plots for the unitary 
simulated pion on the m/ = 0.004, = 0.03 ensemble in figure [I9]and for the kaon in Figure l20l 
The analysis is performed as a simultaneous constrained fit to the five pseudoscalar channels as 
for the 24 3 ensembles (see SectionHII]). The fits are performed between t m [ n =12 and f max = 51. 
We give an example of the reweighted m/ ; dependence of the unitary pion and kaon masses in 
figure I2T1 

Table IXIXI contains the measured pseudoscalar masses and decay constants at the simulated 
strange-quark mass m/, = 0.03. Reweighting to the estimated physical strange-quark mass 
m/, = 0.0273(7), we obtain the masses and decay constants of the pions and kaons in Tables IXXl 
and IXXII respectively. 

We use a gauge fixed box source of size 24 for the Q. baryon using the same configurations as for 
our pion measurements with valence strange-quark masses m x = 0.03 and m x = 0.025 to enable 
an interpolation to the physical strange-quark mass. We display the fit to the m x = 0.03 £1 baryon 
mass on the m/ = 0.004, rrif, = 0.03 ensemble in figure[22l along with the dependence of this mass 
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FIG. 19: Effective pion masses from the PP LW correlator (top left), PP WW correlator (top right), AP LW 
correlator (center left), AP WW (center right) and AA LW correlator (bottom). Note the different vertical 
scale for the WW correlators. The horizontal bands represent the result for the mass from a simultaneous 
fit. 



on the dynamical strange mass under reweighting. We take our fitting range between t m i n = 7 and 

tmax =13. 

The results for the masses of the £1 baryon and the scales r$, r\ and ri/rrj are given in Table lXXlTl 
and IXXHTl respectively, tq and r\ were determined again using Wilson loops formed from prod- 
ucts of temporal gauge links with Coulomb gauge-fixed closures in spatial directions, with an 
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FIG. 20: Effective kaon masses from the PP LW correlator (top left), PP WW correlator (top right), AP LW 
correlator (center left), AP WW (center right) and AA LW correlator (bottom). Note the different vertical 
scale for the WW correlators. The horizontal bands represent the result for the mass from a simultaneous 
fit. 



exponential fit from t = 4 to t = 8 and the resulting potential fit to the Cornell form in the range 
r = 2.45 — 10. An example of the fit to the time dependence of the Wilson loops at the physical 
strange-quark mass is given in Figure [23l This figure also shows the fit to the potential. On these 
ensembles, the strange-quark mass dependence of tq and r\ can be resolved within the statistics, 
but remains small. 
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FIG. 21: We illustrate the m/, dependence of the unitary pion (left panel) and kaon (right panel) masses 
on the mi = 0.004, 32 ensemble. The values are obtained by reweighting around the simulated value 
(m h = 0.03). 
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FIG. 22: We display the fit to the Q. baryon mass with valence strange mass m x = 0.03 on the m/ = 0.004, 
nth = 0.03, 32 3 ensemble showing the quality of the fit with our box source (left panel). We also show the 
weak dependence of the Q. baryon mass with fixed valence mass m x = 0.03 on our simulated my, inferred by 
the reweighting procedure on the mi = 0.004, 32 3 ensemble (right panel). 

A. Nucleon and A Masses 



Baryon effective masses from the 32 3 ensembles are plotted in Fig. [24] and [251 The Gaussian- 
source correlators give good effective-mass signals, while the wall-source correlators are much 
noisier; indeed it is hard to identify a plateau in effective mass signals from the latter. While for 
nucleons effective mass signals from the wall-source seem to eventually settle at the same values 
as from Gaussian source correlators, for the A baryons a plateau cannot be identified from the wall 
source except for the lightest up/down mass. Nevertheless fully correlated fits using two expo- 
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in. 



0.03 0.03 
0.025 0.03 
0.008 0.03 
0.006 0.03 
0.004 0.03 
0.002 0.03 
0.025 0.025 
0.008 0.025 
0.006 0.025 
0.004 0.025 
0.002 0.025 
0.008 0.008 
0.006 0.008 
0.004 0.008 
0.002 0.008 
0.006 0.006 
0.004 0.006 
0.002 0.006 
0.004 0.004 
0.002 0.004 
0.002 0.002 



(0.004) mxy (0.006) (0.008) 



0.3212(3 
0.3073(3 
0.2561(3 
0.2496(3 
0.2430(4 
0.2363(5 
0.2930(3 
0.2392(3 
0.2323(3 
0.2252(4 
0.2180(4 
0.1708(3 
0.1610(3 
0.1506(3 
0.1395(4 
0.1505(3 
0.1393(3 
0.1271(4 
0.1269(4 
0.1133(4 
0.0976(4 



TABLE XIX: Pseudoscalar masses 
simulated strange-quark mass (rrih = 



0.3216(2 
0.3078(2 
0.2565(2 
0.2500(3; 
0.2434(3; 
0.2367(3; 
0.2934(2 
0.2396(2 
0.2327(3; 
0.2256(3; 
0.2184(J 
0.1714(2 
0.1616(3; 
0.1513(3; 
0.1403(3; 
0.1512(3; 
0.1400(J 
0.1280(3; 
0.1278(3; 
0.1144(3; 
0.0989(4 

mxy (mi) and 
= 0.03). 



0.3224(3 
0.3086(3 
0.2579(4 
0.2516(4 
0.2452(5 
0.2388(6 
0.2943(3 
0.2410(4 
0.2342(4 
0.2273(5 
0.2203(5 
0.1727(4 
0.1629(4 
0.1526(4 
0.1417(4 
0.1525(4 
0.1413(4 
0.1293(4 
0.1291(4 
0.1156(4 
0.1001(5 



U (0.004 



0.0801(3 
0.0786(3 
0.0723(3 
0.0715(3 
0.0707(3 
0.0701(3 
0.0770(3 
0.0709(3 
0.0701(3 
0.0693(3 
0.0686(3 
0.0649(3 
0.0641(3 
0.0633(3 
0.0625(3 
0.0633(3 
0.0624(3 
0.0615(3 
0.0614(3 
0.0605(3 
0.0595(3 



Ay (0.006 



0.0804(3 
0.0789(3 
0.0729(3 
0.0721(3 
0.0714(3 
0.0709(4 
0.0775(3 
0.0715(3 
0.0707(3 
0.0700(3 
0.0695(3 
0.0657(3 
0.0648(3 
0.0640(3 
0.0634(3 
0.0640(3 
0.0632(3 
0.0624(3 
0.0623(3 
0.0614(3 
0.0603(3 



/xy(0.008 



0.0809(3 
0.0794(3 
0.0738(3 
0.0731(3 
0.0725(3 
0.0723(4 
0.0780(3 
0.0724(3 
0.0717(3 
0.0711(3 
0.0708(4 
0.0666(3 
0.0659(3 
0.0651(3 
0.0646(4 
0.0651(3 
0.0643(3 
0.0637(4 
0.0634(3 
0.0627(4 
0.0617(4 



the decay constants fxy(mi) on the 32 3 ensembles at the 



nentials to represent the contributions of the ground and first-excited states can be performed for 
both the nucleon and A, yielding the results summarized in Table IXXIVl In addition to this fully- 
correlated two-exponential fit, we have tried two other fit methods: uncorrelated and bootstrap 
correlated with frozen correlation matrix Q40Q . While those earlier analysis were conducted on 
smaller statistics, they agree with the two-state fully correlated fits within two standard deviations 
(see Table IXXVl ) We use the results from the two-state fully correlated fits as our best values of 
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m x m y 


m xy (0.004) 


(0.006) 


(0.008) 


(0.004) 4. (0.006) fxy (0.008) 


0.008 0.008 


0.1706(3) 


0.1711(3) 


0.1725(5) 


0.0645(3) 


0.0653(3) 


0.0662(4) 


0.006 0.008 


0.1608(4) 


0.1613(3) 


0.1628(5) 


0.0636(3) 


0.0645(4) 


0.0654(4) 


0.004 0.008 


0.1503(4) 


0.1510(3) 


0.1526(5) 


0.0628(4) 


0.0636(4) 


0.0647(4) 


0.002 0.008 


0.1392(4) 


0.1401(3) 


0.1417(5) 


0.0620(4) 


0.0630(4) 


0.0641(4) 


0.006 0.006 


0.1503(4) 


0.1509(3) 


0.1524(5) 


0.0628(4) 


0.0636(4) 


0.0646(4) 


0.004 0.006 


0.1390(4) 


0.1398(3) 


0.1414(5) 


0.0619(4) 


0.0628(4) 


0.0638(4) 


0.002 0.006 


0.1268(4) 


0.1278(3) 


0.1295(5) 


0.0611(4) 


0.0620(4) 


0.0632(4) 


0.004 0.004 


0.1267(4) 


0.1276(3) 


0.1292(5) 


0.0609(4) 


0.0618(4) 


0.0630(4) 


0.002 0.004 


0.1131(4) 


0.1142(4) 


0.1158(5) 


0.0601(4) 


0.0610(4) 


0.0622(4) 


0.002 0.002 


0.0974(4) 


0.0988(4) 


0.1003(5) 


0.0590(4) 


0.0598(4) 


0.0612(5) 



TABLE XX: Pion masses m^m/) and decay constants /^(m/) computed on the 32 3 ensembles at the 
physical strange-quark mass m/, = 0.0273(7). 



m x 


m xh (0.004) m xh (0.006) m xh (0.008) 


f xh (0.004) f xh (0.006) fxh (0.008) 


0.008 


0.247(2) 


0.247(3) 


0.249(3) 


0.0712(4) 


0.0718(5) 


0.0727(5) 


0.006 


0.240(2) 


0.240(3) 


0.242(3) 


0.0703(4) 


0.0710(5) 


0.0720(5) 


0.004 


0.233(3) 


0.234(3) 


0.235(3) 


0.0695(4) 


0.0703(5) 


0.0713(5) 


0.002 


0.226(3) 


0.227(3) 


0.229(3) 


0.0687(5) 


0.0698(5) 


0.0710(6) 



TABLE XXI: Kaon masses m x h(nn) and decay constants fxy(mi) on the 32 3 ensembles at the physical 
strange-quark mass m^ = 0.0273(7). 

the baryon masses. They also broadly agree with an independent analysis of baryon masses from 
our ensembles by the LHP collaboration 114 ill within two standard deviations. 

V. COMBINED CONTINUUM AND CHIRAL FITS 

We now turn to the main objective of this paper which is to use the results obtained on the 24 3 
and 32 3 ensembles, as discussed in the previous two sections, to determine physical hadron and 
quark masses and mesonic decay constants in the continuum limit, for physical values of the light 
and strange quark masses. Since we are reporting our first results obtained at a second lattice 
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My 


m h 


m n (0.004) m^(0.006) m n (0.008) 


0.03 


0.03 


0.760(2) 0.765(2) 0.766(3) 


0.025 


0.03 


0.733(2) 0.739(2) 0.740(3) 


0.0273 0.0273 


0.743(6) 0.749(5) 0.753(4) 



TABLE XXII: Omega baryon masses on the 32 3 ensembles at the simulated strange quark mass nth = 0.03 
(first two rows) and at the physical strange-quark mass (third row). 



Quantity 


m h = 0.03 
G(0.004) 2(0.006) G(0.008) 


m h = 0.0273 
2(0.004) Q(0.006) G(0.008) 


ro 


5.52(2) 5.50(2) 5.53(2) 


5.52(2) 5.52(2) 5.55(2) 


n 


3.738(9) 3.718(8) 3.707(9) 


3.754(12) 3.728(9) 3.723(10) 


n/ro 


0.678(2) 0.676(2) 0.670(2) 


0.680(2) 0.675(2) 0.670(2) 



TABLE XXIII: The quantities ro, n and n/ro at the simulated (m/, = 0.03) and physical (m/, = 0.0273) 
strange quark masses on the 32 3 ensembles. Q(mi) denotes the quantity measured with light-quark mass 
mi. 

spacing, we present a careful discussion of our approach to taking the continuum limit and the 
relation between evaluating the continuum limit and determining the physical quark masses. We 
start in Section IV Al with a discussion of what we mean by a scaling trajectory and explain in 
some detail the choice of scaling trajectory which we use in the following. In Section IV Bl we 
describe our power counting scheme, in which we treat the 0(a 2 ) terms in our two ensembles and 
the NLO terms in SU(2) chiral perturbation theory as being of comparable size. In order to gain 
insights into the uncertainties associated with the chiral extrapolation, in addition to SU(2) chiral 
perturbation theory, we introduce an analytic ansatz which is a simple first-order Taylor expansion 
in the light-quark mass. This is explained in Section IV CI We then discuss the specific fitting 
procedure which implements this power counting strategy in Section IVDl and in Section IVEl we 
present and discuss the results. 

A. Defining the scaling trajectory 

Although ultimately we will combine the continuum and chiral extrapolations by performing 
global fits as described in sub section IV A 3 1 and in the following subsections, we start by focussing 
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FIG. 23: The effective potential of the Wilson loops with a spatial extent of r = 2.45 on the ni[ = 0.004 
ensemble at the physical strange-quark mass, overlaid by the fit to the range t = 4 — 8 (left panel). The right 
panel shows the static inter-quark potential V(r) on this ensemble, again at the physical strange-quark mass, 
as a function of the spatial extent of the Wilson loops, overlaid by the fit to the Cornell form over the range 
r = 2.45 - 10. 
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FIG. 24: Nucleon effective mass plots from the 32 3 ensembles. 

on the approach to the continuum limit and discussing the definition and choice of scaling tra- 
jectory. For the purposes of this subsection we imagine that we can perform lattice computations 
for any choice of quark masses and envision performing a series of lattice simulations for a range 
of values of /3, the inverse square of the bare lattice coupling. As /3 — > °° the lattice spacing, 
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FIG. 25: A effective mass plots from the 32 3 ensembles. 

mi N A 

0.004 0.468(6) {4-20} 0.596(15) {4-15} 

0.006 0.498(4) {4-20} 0.615(9) {4-15} 

0.008 0.521(4) {4-20} 0.639(10) {4-15} 

TABLE XXIV: Nucleon and A masses in lattice units from the 32 3 ensembles obtained by two-exponential 
correlated fits to Gaussian-source correlators. {} denotes fit range. 

measured in physical units, will vanish along with all discretization errors. We refer to such a 
one-dimensional path through the space of possible lattice theories as a scaling trajectory. For 
2+1 flavor QCD we must vary the bare lattice mass m M( /(/3) of the up and down quarks and m s (fi) 
of the strange quark so that this trajectory describes physically equivalent theories up to order a 2 
errors. The functions m IK /(/3) and m s (/3) can be determined by requiring two mass ratios (or two 
other dimensionless quantities) to remain fixed as /3 varies. Because of the presence of 0(a 2 ) 
discretization errors, using a different pair of mass ratios will yield a different trajectory of lattice 
theories, whose low-momentum Green's functions will be equivalent to those of the first up to 
0(a 2 ) corrections. 

n 

In ref. Hfl, where we obtained results from simulations at a single value of p, we found that using 
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mi full corr - uncorr. bootstrap a LHP b 
0.004 0.477(4) 0.465(5) 0.469(4) 0.474(4) 
0.006 0.498(2) 0.486(10) 0.489(7) 0.501(2) 
0.008 0.517(3) 0.524(4) 0.5254(16) 0.522(2) 

TABLE XXV: Comparison of nucleon mass results from different analyses on the same 32 3 ensembles. Su 



perscript a denotes Ref . 



40(1 . where a frozen correlation matrix was used and superscript b denotes Ref. 14111 . 



the masses of the % and K mesons and the Q. baryon to determine the lattice spacing a and the bare 
values of m M j and m s was an effective procedure. A natural choice of scaling trajectory would 
therefore be to keep the ratios m K /m.Q_ and m^/m^ fixed as /3 varies. Thus these ratios would 
be chosen to take their continuum values for all /3 with no a 2 corrections. This choice of scaling 
trajectory then fixes the functions w M ^(j8) and m J (/3). In addition, we will identify an inverse 
lattice spacing, expressed in GeV, with each point on this scaling trajectory. To do this we use the 
mass of the Q. baryon and define l/a = 1.672//«^ GeV where 1.672 GeV is the physical mass 
of this baryon and is the mass of the £l~ as measured along our trajectory in lattice units. 
Having defined the scaling trajectory and determined the lattice spacing at each /3 by fixing the 
ratios m^/m^, rnx/m^ and the mass of the Q. baryon to their physical values, we are in a position 
to make predictions for other physical quantities. The results obtained at a particular value of /3 
will differ from the physical ones by terms of 0(a 2 ) . We imagine eliminating these artefacts by 
extrapolating results obtained at several values of /3 to the continuum limit. In order to discuss this 
continuum extrapolation it is convenient to introduce some notation. Let us assume that we have 
performed lattice calculations at a series of /V values of /3, {/3 e }i< e <jv corresponding to points 
along the scaling trajectory defined above (in the present study N = 2). This will determine a 
series of bare quark masses = m^(/3 e ) where / = u d or s. On each of the lattices we compute 
a number of physical quantities, e.g. the kaon leptonic decay constant /£, and our prediction for 
the physical value of fx is the value obtained by extrapolating to the continuum limit. 
Of course, as already mentioned above, the scaling trajectory and the assigned value of the lattice 
spacing at a particular /3 are not unique. Had we used three different physical quantities to calibrate 
the lattice at each /3 and then used the resulting bare quark masses and lattice spacing to compute 
m K /msi, mfc/ina and the mass of the Q. baryon, we would find results which differed from the 
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physical ones by terms of 0(a 2 ). Although there is a choice of the quantities used to define and 
determine the scaling trajectory and the value of the lattice spacing at each /3, for a 2+1 flavor 
theory the number of conditions is always 3iV, where N is the number of different /3 values used in 
the simulations and the factor 3 corresponds to the fact that at each /3 there are three parameters, 
the bare masses m u d and m s and the lattice spacing a. 

In the above presentation we have tried to provide a pedagogical introduction to the determination 
of scaling trajectories and chose to decouple issues related to the extrapolations in the mass of the 
light quark (chiral extrapolations) from the discussion. Of course, in practice at present we are 
unable to perform simulations at physical quark masses, i.e. with masses which give the physical 
values of m K / mci and itik / mci, and so chiral extrapolations are necessary. It will therefore be useful 
in the following to discuss the scaling behavior of a general 2+1 flavor theory in which the masses 
of the pion and kaon differ from those in Nature. Following the conventions defined elsewhere in 
this paper, we will use m/ and for the quark masses in the DWF lattice action which correspond 
to the usual ud and s quarks, and m/ and fhf, for the corresponding multiplicatively renormalizable 
bare quark masses m; = m/ + m KS and = m/, + m res specific to the DWF action. In the next 
subsection we review the origin of the a 2 errors as described by the Symanzik effective theory for 
DWF and in the following subsection present our treatment of scaling for this more general theory. 

1. Symanzik effective theory and a 2 — > extrapolation 

Symanzik's effective theory provides a powerful framework in which to discuss the approach to 
the continuum limit. For any finite value of /3 we expect the low-momentum Green's functions in 
our lattice theory to agree with those in a corresponding effective continuum theory. The effective 
action for this theory contains not only the usual dimension-3 and 4 terms standard in QCD but also 
higher-dimension operators. If the quark masses and the coefficients of these higher-dimension 
operators are properly chosen then the low-energy Green's functions of the lattice and effective 
theories will agree through 0(a d ~ 4 ) provided the effective theory includes all necessary terms of 
dimension up to and including d. This implies that the low-energy Green's functions of the lattice 
theory and the usual continuum theory will differ by the matrix elements of these dimension-5 and 
higher operators which of course are not present in the standard continuum theory. 
For the domain wall fermion calculation presented here the leading corrections come from opera- 
tors of dimension 6. While the dimension-5 Pauli term qo^ v F^ v q is present, its chiral properties 
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imply that it is generated by chirality violation due to propagation between the left and right do- 
main walls. This same residual breaking of chiral symmetry gives rise to the residual mass m res , 
the coefficient of the dimension-3 mass term which remains when the input quark mass is set 
equal to zero. The largest value for m res found in our current calculation, m res = 0.003152(43), 
is suppressed from unity by more than two orders of magnitude. Since a similar suppression for 
this dimension 5 operator is expected, the combination of chiral symmetry and the small value of 
^Aqcd ~ 0.2 suggest this term can be ignored and that the largest finite lattice spacing errors that 
we should expect are 0(a 2 ). 

We require that for our choice of scaling trajectory the matrix elements of these 0(a 2 ) Symanzik 
terms behave as a 2 , allowing a linear extrapolation in a 2 to give the continuum limit. This implies 
that the coefficients of these operators remain reasonably constant along our trajectory. This is 
typically achieved by varying only /3 and quark masses along the trajectory so the only variation 
in the coefficients of these 0(a 2 ) terms comes from the variations in /3 which are quite small in 
present scaling studies [|80D. 



2. Scaling and the quark masses 

In the present calculation we obtain results using a number of light-quark masses, all of which 
are significantly larger than the physical quark masses that were used in the introductory remarks 
above to describe a physical scaling trajectory in which m K /mci, mic/mo, and ma were fixed at 
their physical values. However, we can easily generalize our notion of a scaling trajectory to 
include families of choices for the parameters (/3,m/,m/ 2 ) for which, in an obvious notation, the 
ratios mu/m^h and mih/mhhh are held fixed. In the language used earlier, we require that the N 
triplets of parameters (/3 e , m^, m^), 1 < e < N, lie on the same scaling trajectory if 

ra//(/3 e ,m*,me) 

ra// ? (/3 e ,m*,m<;) 
mhhh(P e ,rn^ml) 

for each pair e and e ; . The ratio of lattice spacings for such a pair would be defined as 

a*_ = m hhh (j5 e ,mlml) ^ 
« e ' m hhh (j3 e ',mf ,m£') 

The scaling trajectory determines two functions m/(/3) and /«/ i (/3), where these bare masses are 
non-trivial functions of /3. While a portion of their /3 dependence should reflect their naive mass 



w//(/3 e ,mf ,m% ) 
m h hh(^ e ',mf,m^) 

mi h (p e \mf ' ,mf) 
m hhh (l3 e \mf ,mi) 



(26) 
(27) 



47 



dimension, these quantities also carry a logarithmic dependence on a characteristic of the anoma- 
lous dimension of the mass operator qq in QCD. Thus, even when expressed as dimensionless 
ratios, e.g. mi(fi)/ma and fhh{fi)/ m Q.i these parameters will have singular continuum limits (in 
fact, the sign of the anomalous dimension of qq is such that these ratios vanish in the continuum 
limit). 

The mass parameters m/ and are short-distance quantities whose definition is free of infrared 
singularities. For example, they could be specified by examining high-momentum, infra-red safe 
Green's functions with no need to compute low-energy masses which are dependent upon the 
low-energy, non-perturbative behavior of QCD. While the individual masses m/(j6) and m/,(j8) do 
not have a continuum limit, both the naive and anomalous scale dependence cancels in their ratio 
m/(/3)/m/ 7 (/3), which is well-defined in the continuum limit and agrees with the corresponding 
ratio in conventional renormalization schemes, such as RI/MOM or MS. 

Let us now assume that we have performed lattice calculations at a series of N values of /3, 
{P e }i<e<N, corresponding to points along the scaling trajectory defined above. This will de- 
termine a series of quark masses = m/-(/3 e ) where f = 1 ox h. It is natural to introduce a series 
of factors which relate the lattice spacings and quark masses between these N ensembles. For con- 
venience, we identify a primary ensemble 1, and introduce 3(N— 1) factors relating each ensemble 
e to the ensemble 1 as follows: 



cf m 



hhh 



1 m 



Z f = £>e for f = ^rh. (30) 

3 K m f 

Since the ratio m//m/ 2 is well-defined in the continuum limit, the corresponding ratio for each 
of these ensembles rhj/m* differs from that limit by a term proportional to (a 6 ) 2 . This 0(a 2 ) 
correction represents the discrepancy between our choice of scaling trajectory with mu/mih fixed 
as we vary /3 and an alternative choice where instead inj/m* is held fixed. Since these trajectories 
differ at 0(a 2 ), we expect that 

^-limf^)(l+c„ 7 (A QCD a e ) 2 ). (31) 
m* 0->oo \m h {p)J 

The term proportional to c m arises from the shifts in m 2 j and mf h caused by the first-order effects of 
dimension-6 terms in the Symanzik effective action. While c m must vanish as in* — > m*, we prefer 
not to write c m as proportional to the difference — because of possible non-analytic terms 



48 



in the quark masses (e.g. possible logarithms of m®) that may appear in the low-energy matrix 
elements of these dimension-6 operators. If we divide Eq. (I3TT) evaluated for our primary ensemble 
1 by the same equation applied to the ensemble e and Taylor expand in the lattice spacing, we 
obtain the following useful relation between Zf and Zf: 



implying the 2 (TV — 1) Z factors associated with the quark masses actually depend on N quantities 
through order a 2 (e.g. we can take the (N — 1) Zf and c m as the independent quantities). The 
constraints implied by Eq. (l32l) do not simplify the N = 2 case addressed in the present paper where 
we would simply be trading the two parameters Z\ and Z ; 2 for the alternative pair of parameters Zf 
and c m . 

Equation (l32l) provides an explicit estimate of how scaling violations revise the standard expecta- 
tion that all quark masses will scale with a common Z factor as the cut-off is varied. As we will see 
from our simulation results presented below, the terms proportional to c m are small and difficult to 
resolve from zero given our statistical errors. 

Since we are now using formulae in which the lattice spacing a e appears alone rather than in a 
ratio, e.g. as a e /a e ' , it may be useful to explain how we intend this is to be determined. It is 
natural to start by considering the physical scaling trajectory discussed in Section lVAl on which 
fnil/mhhh = mn/msi and m//,/m^ = m^/m^. For this physical trajectory, the actual value of 
the Omega mass measured in GeV can be used to define the lattice spacing for any point /3 e on 
that trajectory using a e = m^/( 1.67245(0.29) GeV). In our present study, in order to reach the 
physical trajectory a chiral extrapolation must be performed from the quark masses used in our 
simulation. Ultimately of course, when we present results for dimensionful quantities in physical 
units, it will be necessary to perform the chiral extrapolation and this is the subject of the following 
subsections. For the present discussion of scaling it is sufficient simply to imagine that the lattice 
spacing has been determined in this way and this is the most straightforward way of interpreting 
the 0((a e ) 2 ) terms appearing in equations in this subsection. We stress however, that even this is 
not strictly necessary. We can consider a scaling trajectory defined by fixed, but unphysical, values 
of mu/mhhh and mih/mhhh and define the lattice spacing by assigning an arbitrary value to Mhhh> 
the mass of the hhh baryon on the trajectory in "physical" units, a e = mf lhh / 'M^hh- While the value 
of a e defined in this way depends, of course, on the choice of M/,^, this arbitrariness is simply 
absorbed by a change in constants such as c m in (|3Tb . For the discussion in this subsection it is 




(32) 
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sufficient to note that such a definition of the lattice spacing is possible in principle, the numerical 
determination of a e does not actually have to be performed. 

In the analysis to follow we will examine a family of nearby scaling trajectories in which m/ 
and vary over limited ranges (specifically, m/ varies up to about 0.013 on our coarser lattice 
and m/j varies by up to 20% around m s ). Consider two such trajectories, defined by keeping the 
ratios mu/mhhh an d mih/ m hhh fixed along each trajectory, but taking different values on the two 
trajectories. Let mu/m^h = r u and m lh /m hhh = r {h on the first trajectory and mn/m hhh = r' u and 
m lh/ m hhh = r 'ih on the second. As /3 — > °°, the ratio of bare quark masses on the two trajectories 
will approach a limit up to 0(a 2 ) corrections: 

where f=l or h, and m*(r//,r//j) and m^(r/;,r//j) and ^( r /Z» r /ft)) are the values of the 

bare quark masses on ensemble e such that mu/mh.hh = r u an d mih/ m hhh = r ih ( m u/ m hhh = r \\ 
and m lh /m hhh = r' lh ). The ratios R a = m\ hh {m^(ri h r lh ),m\{ri h ri h )) /m^(mf (r;/,^),^^//,^)) 
and < = m^ A (mf(r; z ,4),mJ(^,r;j)/m^ M (m^^,4),m^r z ' z y tt )) each describe the change in 
lattice scale as the bare coupling changes from /3 1 to /3 e . In the limit of small bare coupling, this 
change of scale can be determined entirely from the short-distance part of the theory and must be 
the same for our two trajectories up to order a 2 corrections since these two trajectories differ only 
in the choice of quark masses. Thus we can write 



a T = l+d a A 2 QCD ((a*) 2 -(a 1 ) 2 ) (34) 



where we have explicitly represented the fact that each ratio and hence the ratio of ratios must 
approach unity as a e — > a 1 . Both the coefficients d m j and d a will vanish when the primed and 
unprimed trajectories that are being compared become identical. 

Taking the ratio of two versions of Eq. (|33~1) . one for /3 e and the other for our primary ensemble 
^and using Eq. ([341) , we obtain an expression for the change in the factors Zf between these two 
trajectories: 

( K \ + (d mJ + d a )A 2 QCD [(a 1 ) 2 -(a e ) 2 ]y (35) 

Since the changes in fh[ and m h between these two trajectories which we wish to compare are small, 
the resulting coefficients d m j and d a will also be small and we will neglect the 0(a 2 ) correction 
on the right-hand side of Eq. (1351) . Thus, we will use the same values for Z/ and Z/ 7 for this family 
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of nearby trajectories, i.e. we drop lattice artefacts proportional to m } and (m h — m s ) and so neglect 
the mass dependence of Z/ and Z/, in this limited range of masses. In the following we will refer 
to this range for mi and as their "allowed range". 



3. Fitting strategies 

We exploit the above relations between numerical results obtained at the two values of /3 for which 
we have performed simulations in two ways. The first we label the "fixed-trajectory" method. In 
this approach we determine R a , Z\ and Z/, by matching results obtained at a single pair of equivalent 



quark masses 118 111 . For example, the masses used at one value of /3 may correspond to values at 
which a simulation was actually performed. The corresponding set of masses for the other /3 might 
be determined by linear interpolation to make the two ratios mn/m^ and mi^/m^ agree with 
those on the first ensemble. The ratio of lattice spacings and the two Zf factors are then determined 
from Eqs. (l29l) and d30l) . It will be important to recall that Z/ and Z/ 3 are constant in the allowed 
range of quark masses. Finally, knowing the three factors R a , Zi and Z h we make a common fit to 
the mass dependence of physical quantities computed for both values of /3. 
In the final step, we adopt an ansatz for the mass dependence that is expected to be accurate 
both for the points in our calculation and for the physical values to which we wish to extrapolate, 
specifically a NLO chiral expansion about the chiral limit or a simple Taylor expansion about 
the physical point. Each ansatz for the continuum theory, when combined with the three scaling 
factors R a , Z\ and Z/ 7 and with any required a 2 corrections, will then provide a set of formulae 
which should describe all of our data for both /3 values. For example, in the chiral fits described in 
the next section we can use a common set of Low Energy Constants (LECs) to fit both sets of data 
provided we scale the values used on one set by the required factors of R a , Z/ and Z/, before we 
use them on the other. Where explicit 0(a 2 ) terms are required, these can be added with unknown 
coefficients which are also scaled appropriately between our two values of /3 . In such a combined 
chiral and a 2 expansion we adopt a power counting scheme, described below, so that only effects 
of a similar minimum size are consistently included. 

During the initial process of determining R a , Z\ and Z/, we cannot assign a physical value to the 
lattice spacing. The original trajectory being used does not correspond to physical masses so no 
notion of "GeV" exists for that case. Of course, the further fitting to the quark mass dependence of 
the two ensembles is introduced to allow extrapolation to physical values for the ratios mn/m^ 
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and mih/mhhh- When is evaluated at this same physical point, its value can be compared with 
1.672 GeV to determine the lattice scale. 

This fixed trajectory method is intended to cover a wider range of possible scaling trajectories than 
the example discussed above where the trajectory passes precisely through one of the simulation 
points. If we wish, we can adopt an ansatz for the quark mass dependence of m n , yyik and and 
perform this fixed trajectory scaling with the parameters R a , Z/ and Z/, allowed to vary and fix their 
values from Eqs. (l29l) and (|3Q|) at values of m/ and m/ 7 for which the ratios mn/mhhh and mih/mhhh 
take their physical values. 

The second approach, termed "generic scaling", introduces the factors R a , Z/ and Z/ 7 as parameters 
into the ansatz being used to fit the quark mass dependence. In this approach we perform a fit to all 
our data for m K , mx and over a range of quark masses for which the fitting ansatz is accurate 
and for which the use of fixed values for R a , Zi and Zj 7 is legitimate. In this generic scaling ap- 
proach, our choice of scaling trajectory with fixed hadron mass ratios mu/mhhh and mih/mhhh an d 
with nikhh determining the lattice scale is realized somewhat indirectly. The three conditions asso- 
ciated with this choice of scaling trajectory are realized by omitting possible a 2 corrections from 
the expressions used to fit mu, mih and mhhh- The resulting trajectory can therefore be interpreted 
as being the one along which the masses of the pion, kaon and Q-baryon take their physical values, 
as was the case in the discussion of Section IVAl The difference of course, is that whereas in Sec- 
tion IV Al we envisaged (unrealistically at present) being able to simulate directly at the physical 
value of mi, we now reach the physical point after an extrapolation in quark masses. The detailed 
discussion of the ChPT functions used in describing the quark mass dependence of the pion and 
kaon masses is given in Sub section IV Bl and those for the analytic ansatz in Sub section |VC] below. 
However, both our ChPT and Taylor expansion ansatze stipulate that to the order being studied 
mhhh is a linear function of rh\ and 7%. It is instructive to explore this case here. 
Included among the equations used to determine the low energy constants and the scaling factors 
R a , Zi and Zh are two equations for mhhh on our two ensembles: 

m\ hh {m h m h ) = ml hh (0,m h0 ) +<4 nm; m/ + cJ, Q m*(^A -m h0 ) (36) 



—mlhh(RaZimi,RaZ h mh) 



(37) 



Ra 

= J- m hhh(°> ™ho) + c} nnmi R a Zim, + c} namh [R a Z h mh - mho) 

Here 1 is our primary ensemble, for us that is the one with — 2.25 and the 32 3 x 64 volume, while 
the second ensemble is the one with the coarser lattice spacing and is labeled 2. m* W! (m/,m/,) are 
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the /z/z/z-baryon masses corresponding to bare-quark masses m/ and on ensemble e. Although 
we have written m^o as a general constant, we have in mind to use the equations with m/ ? o in the 
allowed range of the physical bare strange quark mass in the primary ensemble. Equations (l36l) 
and (|37T) define the three constants m^(0,m/jo), cj, m and c} nQpih which are related to the physical 
Q. mass and its "physical" dependence on the quark masses. The absence of 0(a 2 ) corrections 
to Eqs. (l36l) and (1371) implements our choice that is being used to set the scale and hence by 
construction contains no finite lattice- spacing errors. While part of a larger set of equations which 
are being used to determine the low energy constants as well as R a , Z\ and Z h , the leading order 
effect of these two equations is to determine R a . Note that this is identical to imposing Eq. d29l) 
in the fixed trajectory method at the point mi = 0, m/ 7 = mj 7 Q. Since the variation of R a as m/ and 
m/, change over their allowed range is of the same size as the variation of Z/ and Z/, over this same 
range it can also be neglected, so any particular choice of m/, is equivalent to any other within this 
allowed range. 

The fixed trajectory and generic scaling methods are similar in nature. Both require that an ansatz 
be adopted to allow the quark mass dependence of lattice quantities to be described in order to 
define the scaling parameters R a , Zi and Z/ z and to extrapolate to the physical point. Both assume 
that the scaling relations between the two ensembles defined by R a , Z/ and Z/, hold over the allowed 
range of masses. The fixed trajectory method corresponds most closely to our original definition 
of a scaling trajectory and decouples the matching of the two lattices from the chiral extrapolation. 
It requires however, the introduction of a convenient but arbitrary point at which the matching 
between the two ensembles is performed. The generic method avoids this arbitrary choice and 
applies these assumptions uniformly over the entire range of allowed masses. The fixed trajectory 
method determines R a , Z/ and Z/, in an iterative fashion as explained in Section lVDl The generic 
approach determines the coefficients in the adopted ansatz from a single % 2 minimization. The 
physical quark masses are then determined by inverting the resulting equations which give m n , mx 
and mo, in terms of m/ and m^. 

The detailed discussion and results presented in this paper correspond to the fixed trajectory 
method; fits using the generic scaling approach were performed to monitor the consistency of 
the results and estimated errors. 
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B. Scaling and chiral perturbation theory 

At the start of section lVAl we discussed the continuum extrapolation in an idealized situation in 
which we can perform simulations at any value of the quark mass m/. In reality this is not the case; 
for example, the lightest unitary pion appearing in the current study has mass 290 MeV. In order to 
compare our results with Nature we therefore need to extrapolate to lighter quark masses and this 
was already acknowledged when discussing the fitting strategies in section lV A 3 1 above. We now 
explain how we combine the continuum and chiral extrapolations in global fits. We start in this 
section by using SU(2) chiral perturbation theory for the mass dependence, with the expectation 
that the extrapolation will be made more precise if constrained by the theoretically known behavior 
of QCD in the chiral limit [1]. However, in order to estimate possible systematic errors associated 
with this extrapolation and to obtain a more complete understanding of the implications of our 
calculation, we also examine a simpler analytic extrapolation to physical quark masses [42] and 
this is explained in the following subsection. Although later we will perform extrapolations using 
partially quenched ensembles, for the purposes of this introduction we restrict the discussion to 
the unitary theory in which the valence and sea quark masses are equal. 

We now explain the power counting scheme we employ to identify NLO corrections to the chiral 
and continuum limits. Since the pion mass and decay constant are central to SU(2) ChPT, we 
begin by considering the predictions of continuum NLO ChPT for these two quantities: 

mfi = Xi +Xi ■ < ^( (2IT -4 2) ) +2(2zf -^))xi + J^TpXi^ > ™ 




fu = f + f-{^X'+Lf) Xl -^\og^ \. (39) 

Here mu and //; are the mass and decay constant of the pseudoscalar meson composed of two 
light quarks, /, L4, L5, L(, and L% are the conventional low energy constants and A^ is the usual 
chiral scale. The quantity Xi comes directly from the lowest order chiral symmetry breaking term 
in the effective chiral theory and is proportional to the QCD light quark mass. It is conventionally 
written Xi = 25m/, where B is another low-energy constant. 

We now discuss how we apply these formulae to describe the low energy behavior of lattice the- 
ories which lie on a scaling trajectory. For a sequence of ensembles {e}i< e <# lying on such a 
scaling trajectory not only will the quark masses and lattice units, (m®,m®,a e ) be related, but also, 
when expressed in physical units, the quantities /, L4, L5, and Ls should take the same val- 
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ues up to 0(a 2 ) corrections. The same is true for the renormalization independent combination 
Xi = 1Bm\ (see the discussion below). As detailed in Ref. yj], chiral perturbation theory at finite 
lattice spacing for domain wall fermions involves a simultaneous expansion in the explicit bare 
quark mass, mi, the squared lattice spacing, a 2 , and the residual chiral symmetry breaking arising 
from the finite separation, L s , between the two four-dimensional walls in the fifth dimension. We 
will denote this last quantity by e~^ Ls , suggesting the exponential decrease in such residual chiral 
symmetry breaking found in perturbation theory for DWF. (The actual behavior is a sum of ex- 
ponential and inverse power dependence on L s .) No new terms need to be added to the resulting 
effective low energy theory to describe the resulting Green's functions to NLO in the parameters 
m/, a 2 and e~^ Ls . Thus, we can use equations with the form of Eqs. (1381) and (l39l) to describe 
the lattice results for m// and fu along a scaling trajectory, provided we work to NLO in a power 
counting scheme which treats the quantities Xl/i^nf) 2 , a 2 AQ CD and e~^ Ls as equivalent and keep 
a single power of any of these quantities as a correction. We must now determine how the param- 
eters appearing in these equations must be adjusted to describe lattice results at finite a 2 . 
Since the scale can be freely varied if the other analytic terms are appropriately changed, we 
will choose this quantity to be constant if measured in physical units. Thus, for each point on 
our physical scaling trajectory we will choose A^ = biq • 1/1.672, giving it the value of 1 GeV. 
Because of their proportionality to the NLO factor Xi all of the parameters which appear in the 
large curly brackets on the right hand side of Eqs. (1381) and (l39l) can be given their continuum 
values, dropping possible 0(a 2 ) terms as being of NNLO in our power counting scheme. Thus, 
within those brackets the quantities /, L4, L5, Lg and L%, when expressed in physical units, can be 
given identical values for the ensembles on the scaling trajectory. 

In contrast, when Eq. (l39l) is used to describe our finite lattice spacing results, the LO quantity f e 
determined on ensemble e, expressed in physical units, depends on /3 e . However, it approaches its 
continuum limit with 0(a 2 ) corrections and so we write f e = f + Cf(a e ) 2 . 

Given the definition of a scaling trajectory, the variation of the quantity xf needed to apply Eq. (1381) 
to the ensemble e is actually trivial. Because our choice of quark mass gives the same value for 
inn f° r eac h ensemble e on our scaling trajectory, all of the quantities in Eq. (1381) with the possible 
exception of the xf which we are now considering, are the same when expressed in physical units 
for all points on the scaling trajectory. Thus, xf = 2B e mJ/ (a e ) 2 must be a constant as well, where 
B e and are explicitly left in lattice units. Since we know how the quantities m/ and a 2 are related 
between an ensemble e and our primary ensemble 1, we can determine the N — I constants B e in 
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terms of the single constant B 1 : 

Z e 

B e = -Lb 1 (40) 

R e 

a 

without any a 2 corrections. Because of the complex scaling behavior of the mass, we will treat 
B 1 as one of the LEC's to be determined in our fitting and not relate it to a "physical" continuum 
quantity whose definition would require introducing a continuum mass renormalization scheme. 
We conclude that our lattice results for light pseudoscalar masses and decay constants obtained 
from a series of ensembles {e} can be described through NLO by the formulae: 



ml,) 2 = Xf + Xf { ^^-Lf) + 2(2Lr-^)Jzf + T ^Zflog A 2 
ft = f [l +c f (a'f] + f ■ | £(21® + Lf) X t - g^logf 




with 
,e 



Zf B 1 ^ 



xf = ibr^ («) 



where all quantities in Eqs. (l4TT) and (l42l) are expressed in physical units (except for B 1 and in 
Eq. (1431) which are given in lattice units). 

Two important refinements should be mentioned. First, for the case of a physical scaling trajectory, 
i.e. one which terminates in the physical masses m K , mx and m^, these physical units are naturally 
GeV. However, for other scaling trajectories appropriate "physical" units to use can be those in 
which the Omega mass is unity. Second, for simplicity in Eqs. (|38|) . (|39| ), (14TI) and (l42l) we have 
treated the heavy quark mass as fixed and not displayed the dependence of the quantities /, B, 
L4, L5, Lg and Lg on m/,. In practice we can easily generalize these equations to describe the 
dependence of m// and /// on m/, as well. Provided we limit the variation of to a small range 
about an expansion point fhhQ, this variation can be described by including a linear term in m/, — fhho 
and treating this term as NLO in our power counting scheme. Thus, such extra linear terms will 
only be introduced into the leading order terms in Eqs. (|4TI) and (l42l) . 

Next we present the corresponding formulae for the quantities m# and itiq which are used in the 
determination of the scaling trajectory and in the assignment of a lattice spacing at each value of 

W) 2 = (» W ) 2 +(» W ) 2 {^2f} (44) 
m hhh = m {a) +m {a) c ma , mi xf- (45) 
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Here mS K ^ and mS°^ are the mass of the Ih meson and the hhh baryon respectively in the SU(2) 
chiral limit, i.e. with m/ = 0, for the value of m/, used in the simulation. Similarly the LECs X\ 2 
and Cm^m; depend on m/ z and we are using the notation for the LECs A12 which we introduced 
in Jl]. (Note that c mQmi , whose value is given in Table IXXVIII below, should be distinguished 
from the related parameter cjj, OT which appears in Equations (1361) and (1371) above.) The absence 
of any corrections of 0(a 2 ) on the right-hand sides of Eqs. (|44|) and (|45T) follows from the same 
argument which justified omitting an 0(a 2 ) correction from the right hand side of Eq. (|4Tj) . For 
masses mj and nf h lying on a scaling trajectory the left hand sides of these equations must all be 
the same because of our definition of scaling trajectory. Because of our power counting scheme, 
no a 2 corrections need to be included in the NLO terms proportional to xf on the right hand side 
of these two equations. Therefore the leading order terms and m^> must also be the same for 
all ensembles when expressed in physical units and no 0(a 2 ) correction can appear. As discussed 
above, these equations can be generalized to describe the NLO dependence on m/, varying about 
an expansion point m^. In fact, for the Q. baryon this more general case for Eq. (1451) was described 
in the previous subsection in the equivalent Eqs. (l36l) and (1371) . 

Note that the coefficient of the chiral logarithm in Eq. (|4TI) includes a factor which depends on /, 
the pion decay constant in the SU(2) chiral limit (all other factors of / in Eqs.(l4TT) and (1441) can be 
absorbed into a redefinition of LECs which in any case are determined by fitting). This low energy 
constant / can be determined from the measured values of //; using Eq. (1421) . but to NLO it can 
also be replaced by the measured values of fy. 

As described in Subsection IV A 31 these ChPT formulae can now be used to determine physical 
results in the continuum limit from those obtained on our two lattice spacings. We can employ the 
fixed trajectory method, finding the ratios Z/ and Z/, which relate a specific choice of quark masses 
on one ensemble to those on the other which lie on the same scaling trajectory. The corresponding 
ratio of values of mhhh determines R a . These three quantities then allow a single set of LECs to 
be used to extrapolate the results of both ensembles to the continuum limit and to the physical 
value of the light quark mass using Eqs. (14TT) . (l42l) . (|44l) and (1451) . As a result we learn the physical 
values of m, ( ^(/3 e ), m s (l5 e ) and a e on our two ensembles. In other words, we determine the quark 
masses and lattice spacings for our two ensembles which lie on the physical scaling trajectory. 
Alternatively, we can use the generic fitting approach and introduce the three parameters 
(ZijZfaRa) into the four equations Eqs. (14TT) . (l42l) . (|44|) and (l45l) and obtain a fit to the lattice data 
from both ensembles for which the quark masses lie in the allowed range. The resulting values of 
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the LECs and (Z/,Z/ 17 i? a ) then determine the functions mf,(m/,m^), mf> (m/,m/j) and ra?>, (m/,m/j). 
The physical quark masses on each ensemble, m* rf = m Md (/3 e ) and = m 5 (/3 e ), are then obtained 
by solving the equations: 

mf z (m^,m£) _ ^ ^ mj^g^mg _ m K ^ 



m hhh«di m 1) m n m hhh( m ld' me s) m £i 

where on the right-hand sides the ratios take their physical values. 

Having determined m [K /(/3 e ), m s (/3 e ) and a e as described above, we are in a position to compute 
other physical quantities. For example, at NLO in our power counting the behaviour of the kaon 
decay constant f K is 



f',H = /« fi+c fM K) 2 l + H - (^l^f ' • (47) 



where is the result in the SU(2) chiral limit (m; = 0), ^3,4 are m/, -dependent low-energy 
constants and c^ K ) is a constant. For each /3 e , having determined m s (/3 e ) we measure ff h for 
= m. s (/3 e ) as a function of m/; fit the measured values at all /3 e to determine the LECs and 
Cj(k) in Eq. (PTTT) and finally obtain the physical value of fx by setting a = and m/ = m [(rf . Such a 
procedure is then generalized to the other physical quantities we wish to compute. 

C. Scaling combined with an analytic ansatz for the chiral dependence 

While we know that the ansatz based on chiral perturbation theory described in the previous sub- 
section is valid in the limit of small u and d quark masses, we do not know the precision with 
which it holds over the range of masses which we analyze in this paper (corresponding to data in 
the range 240 MeV < m K < 420 MeV). Indeed it is precise lattice simulations which will answer 
such questions. In order to obtain some understanding of the corresponding systematic uncertain- 
ties, in addition to the procedures based on chiral perturbation theory described in section lVBl we 
consider an ansatz based on a first-order Taylor expansion about a non-zero quark mass, in the 



style of ref. II42L 14311 . Within this approach, since we do not include chiral logarithms, we are not 
able to take the chiral limit and only assume the validity of the analytic ansatz between the physical 
point (to which we extrapolate) and the region where we have data. In this work we only consider 
linear, first-order fits and are therefore insensitive to the choice of expansion point which we take 
to be the same as that at which we match the ensembles when using the fixed trajectory method. 
This simplifies the discussion below of the simultaneous expansion in a 2 and mass differences. 
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Beyond first order, convergence may be improved by considering an expansion point between the 
region in which we have data and the physical point, but this is beyond the scope of our current 
analysis. 

Using the analytic ansatz for m 2 K as a function of the quark mass m q , we find numerically that the 
constant (mass independent) term is consistent with zero, indicating that the tangent of m\(m q ) 
in the unitary case does pass through the origin. Thus, at our statistical precision, no significant 
chiral curvature is needed to satisfy Goldstone's theorem, however we retain the view that we are 
indeed using a model which is valid only in a restricted region of non-zero quark masses. 
Goldstone's theorem also applies in the partially quenched theory and the pion mass vanishes as 
the valence-quark masses are taken to zero while keeping the sea-quark masses fixed. In this case 
however, our linear fit extrapolates to a non-zero pion mass for massless valence quarks, and this 
naturally implies that some form of curvature is required at smaller masses. This is consistent with 
enhanced chiral logarithms in the partially quenched theory. However, the fits do not necessarily 
imply that chiral logarithms at NLO correctly represent the quark-mass dependence between the 
simulated range of masses and the physical point. Instead, in this approach the sum over multiple 
orders of chiral perturbation theory is assumed to be approximated by a linear dependence in the 
relevant range of masses. It is also possible of course that the simulated range of masses is outside 
the useful domain of chiral perturbation theory and that, for example, phenomenological models 
based on combining NLO chiral perturbation theory with arbitrary analytic subsets of terms which 
appear at NNLO and NNNLO are less well motivated than our linear ansatz. 
For m\ and f K it is convenient to define the average valence quark mass m v = m " ^ my . As in 
section lVBl we apply a power counting rule in a double expansion in m x — m m , m y — m m , mi — m m 
and a 2 , where m m is the mass at which we match the ensembles which we also choose to be the 
point around which we perform the Taylor expansion and we recall that m x>y and m/ are the valence 
and sea light-quark masses respectively (here we allow for partial quenching). For the pion mass 
we use the ansatz 

m% = C%* + C'^ (m v - m m ) + C™* (m t - m m ) , (48) 

where we use our standard notation in which the subscripts xy imply that the two valence quarks 
have mass m x and m y respectively. By the definition of our scaling trajectory, there is no 0(a 2 ) 
term at the match point and so there is no correction to C™ 71 . Within our power counting we could 
equivalently use 

m 2 ^ = + C™"m v + (%*mi , (49) 
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where for convenience we redefine C™ 71 between equations (l48l) and ( |49l ). 

In searching for evidence of chiral logarithms it is conventional to plot the ratio m 2 y /m v as a 



function of the quark masses. With the ansatz proposed in Eq. (|49l) 

.2 r m % n m % ■ 



in 



+ C m„ + Z A^L^ (50) 



m v m v m 



V 



and we note that an observed deviation of the mass dependence of from a constant in the 
finite range of quark masses which can be simulated, is not in itself unambiguous evidence of a 
non-analytic structure. 

For decay constants, which do not vanish in the chiral limit, the 0(a 2 ) term are not sensitive to the 
choice of expansion point: 

/xy = C^[l+C f y]+c{^m v -m m )+C^(ffi 1 -m m ) (51) 
= Cfr[l+C f a 2 ]+C{ n m v + C( n mi, (52) 

where again we have redefined Cq between the first and second lines. 

Following a similar argument, at a fixed strange-quark mass, we take the light-quark mass depen- 
dence of the kaon mass and decay constant and the mass of the ^-baryon to be given by 

m xh {a,mi) = C™ K + C™ K m x + C™ K mi , (53) 
f xh (a, mi) = C f *[\ + C fK a 2 ] + c{ K m x + C f 2 K m t . (54) 
m hhh (a,mi) = C^ Q + C* Q m/. (55) 

We stress that the constants C™ n , c[ n , Cf, C™ K , C„ K , Cf K and C'" n implicitly depend on the strange 
quark mass. 



D. Procedure for combined scaling and chiral fitting 

Having introduced the theoretical framework behind our combined scaling and chiral fits in Sec- 
tions |VB] and |VC] we now explain its practical implementation. The formulae given above which 
describe the combined behaviour are valid only for a fixed strange-quark mass and we are faced 
with the problem that the physical strange mass is not known a priori but is an output of the 
calculation. The procedure for performing the combined chiral-continuum fits is therefore neces- 
sarily iterative. As explained in more detail below, we start with some initial values for the lattice 
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spacings and quark masses, perform the fits and then use linear interpolations in m/ 7 to obtain 
updated estimates. The process terminates when the updated estimates converge. During this it- 
erative procedure we use reweighting (see section llTDl) to adjust all pionic observables to the new 
strange-quark mass on each ensemble. For kaon and Q. observables a linear interpolation between 
the unreweighted unitary measurement, and measurements with a second valence strange quark 
(reweighted-to-be-unitary) suffice to obtain that observable for m y = m h = mf uess . 
For the remainder of this subsection we explain further the procedure which we use to match 
lattices with different /3 and present results for the ratios and Z* defined in Eqs. (|29l) and (1301) 
for our ensembles using the fixed trajectory method explained in Section lV A 31 We start by taking 
a specific value of (m/,m/j) M on the ensemble M to which the other ensembles are matched. We 
refer to this as the matching point. The ensemble set M may be the same as the primary ensemble 
1, but does not need to be. As discussed in section IVAl the matching to other ensembles e ^ M is 
performed by requiring that the ratios of hadronic masses and are the same on all lattices 

m hhh m hhh 

at the matching point. Although the final physical predictions do not depend upon the choice of 
matching point, certain choices are favoured due to the quality of the data at the matching point 
and the range over which the data must be interpolated/extrapolated on the other ensembles to 
perform the matching. The ideal point has as small a statistical error as possible and lies within 
the range of simulated data on all of the matched ensembles such that only a small interpolation is 
required. In practice, the errors on the mass ratios at the matching point can be reduced by fitting 
to all partially quenched simulated data on the ensemble set M and interpolating to the matching 
point along the unitary curve. We use linear fitting functions for the light-quark mass dependence 
of the pseudoscalar mesons and the £1 baryon in these short interpolations: 

mly = co + cimi + c v (m x +m y ), (56) 
m xh = do + dimi+d v m x , (57) 
mhhh = e + eim h (58) 

where as elsewhere x,y (I) represent the light valence (sea) quarks and h represents the heavy 
quark. Equations (1561 - (1581 ) are written in lattice units. Although the linear behaviour in Eqs. (|56l ) - 
(1581) is similar to that used in the analytic ansatz, Eqs. (l49l). (1531) and (1551) . we stress that the 
meaning is different. When using the analytic ansatz we assume its validity in the full range 
of masses between the physical ones and those we simulate. Eqs. (l56l) - (1581) on the other hand, 
are only assumed to represent the mass behaviour in the short intervals between the matching 
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and simulated points on ensembles e ^ M, independently of whether we subsequently use chiral 
perturbation theory or the analytic ansatz to perform the chiral extrapolation. 
Once a matching point has been chosen, the matching proceeds as follows: 

1 . For each set of ensembles e ^ M, we perform an independent partially-quenched linear fit 
to the simulated pion, kaon and Omega masses using the forms given in Eqs. (l56l) - (l58l) . 

2. We make a first estimate of the pair of quark masses (m/,m^) e on each ensemble set e ^ M 
that corresponds to the matching point. 

3. We then interpolate the three hadronic masses to the estimated for each value of the 
simulated unitary heavy quark mass. 

4. We linearly interpolate each quantity to the estimated value of mjr. 

5. Next we calculate the ratios R] = and R\ = 

m hhh m hhh 

6. Using the measured slopes of and m^ hh with respect to m®, by comparing R* to the 
corresponding value Rf 1 at the matching point we obtain an updated estimate of m*. 

7. Similarly, by comparing the ratio R* to R™ we obtain an updated estimate of m e h . 

8. With these updated estimates of the quark masses (m/,m/j) e , we return to step[3]and iterate 
the steps until the process converges. 

Once this procedure has converged, we have a set of bare quark masses (m/,m/ 7 ) e which, in phys- 
ical units, are equivalent to the masses {mi,mh) M . Following the discussion in Sec. IV A2l we 
choose a primary ensemble 1 and determine the ratios of quark masses Zj in ensembles 1 and e as 
in Eq. (l30l) with the corresponding ratios of lattice spacing R a given in Eq. (l29l) . 
In the above we assumed that for each ensemble e we had performed simulations at several val- 
ues of m*. In our present study the simulations were performed at a single value of and the 
dependence on the heavy-quark mass is obtained by reweighting as explained in Section llTDl 
The above discussion was deliberately presented in a general case where there are an arbitrary 
number of ensembles. In our case we only have two sets, i.e. the 24 3 and 32 3 lattices. For the 
primary ensemble we choose the finer 32 3 lattice. As we have only one other ensemble set (24 3 ), 
from now on we drop the superscript on the ratios of lattice spacings (R a ) and quark masses (Z/ 
and Z h ). 
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M 


(a mi ) M 


(am h ) M 


(ow/) e (amh) e Z[ 


Z h 




32 3 


0.004 


0.03 


0.00313(13) 0.03812(80) 0.980(15) 0.976(11) 0.7617(72) 


32 3 


0.006 


0.03 


0.00583(12) 0.03839(51) 0.981(9) 


0.974(7) 


0.7583(46) 


32 3 


0.008 


0.03 


0.00860(19) 0.03869(64) 0.979(10) 


0.972(8) 


0.7545(58) 


24 3 


0.005 


0.04 


0.00545(11) 0.03148(51) 0.985(12) 


0.978(9) 


0.7620(57) 


24 3 


0.01 


0.04 


0.00897(18) 0.03074(57) 0.974(11) 


0.968(9) 


0.7517(70) 



TABLE XXVI: Values of the quark mass ratios Z/ and Z/ t and the lattice spacing ratio R a determined by 
matching at five points over both ensemble sets. The quark masses here are quoted without the additive m res 
correction. The ensemble e 7^ M. 
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FIG. 26: Ratios of dimensionless combinations of lattice quantities Q (listed in the figure) between the 32 3 
and 24 3 lattices at the matching point corresponding to ra/ = 0.006, ra/, = 0.03 on the 32 3 lattice. A value of 
unity indicates perfect scaling. The ratios mu/m^h and mujm^h (and consequently ra///ra//,) are defined 
to scale perfectly at these quark masses as a consequence of our choice of scaling trajectory. 



In Table IXXVII we give results for Z/, Z/, and R a obtained by matching at several matching points 
on both ensemble sets M G {24 3 ,32 3 }. Since we prefer to have a matching point within the range 
of simulated data on both ensembles, we can discard the first and last entries in the table. From 
the remaining 3 possibilities, we choose as our final values Z/ = 0.981(9), Z/, = 0.974(7) and 
R a = 0.7583(46) from the second entry with M = 32 3 and (m h m h ) 323 = (0.006,0.03). 
Having chosen to perform the matching of the lattices at the two lattice spacings by requiring 
that mu/mhhh and mih/nihhh take the same values at the matching point, we expect to see lattice 
artefacts in ratios of other physical quantities. This is illustrated in Figure|26]in which we show the 
ratios of several other dimensionless combinations of lattice quantities between the two lattices at 
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the quark masses used in the matching procedure above. The figure shows that we can expect only 
small scaling violations on the order of 1-2% for the other quantities used in our global fits, and 
also confirms that other dimensionless combinations of lattice quantities would be equally suitable 
choices for the definition of the scaling trajectory. 



E. Results of combined scaling and chiral fits 

Using the matching factors Z/, Z/ 2 and R a determined as described in the previous section we are 
ready to perform a simultaneous fit of all our pion, kaon and Q. mass and decay constant data 
to either the NLO forms in chiral perturbation theory, Eq. (l4~TT) to Eq. (l45l) . or the analytic forms 
Eq. ( |49l ) to Eq. (1551) . We also correct for finite volume effects in NLO PQChPT by substituting the 



chiral logarithms with the corresponding finite-volume sum of Bessel functions 114411 . The iterative 
procedure is the same for each of these three fit ansatze. For each iteration i, we: 

1 estimate the physical strange-quark masses, m l s , from the (i — l)th iteration; 

2 interpolate and reweight the data to m\\ 

3 fit the ra x , ra v , ra/ dependence of the light pseudoscalar mass and decay constant; 

4 fit the m x , mi dependence of kaon quantities at m/, —m l s ; 

5 fit the ra/ dependence of the Omega mass for m h — m' s ; 

6 by comparing to the physical values of m K /ma and ra^/ra^, determine the iterated predic- 
tions for the physical strange quark masses m l s +1 



This process is repeated until it converges and a self consistent set of quark masses, lattice spacings 
and results in the continuum limit are obtained. 

For the fits based on NLO chiral perturbation theory we use Eqs. (l4Tj) and (142)) for the pion mass 
and decay constant respectively, and Eqs. (l44l) and (l47l) for the kaon mass and decay constant. 
In our earlier work yj] we found that we had to apply cuts to keep the pion mass below around 
420 MeV in order for NLO SU(2) ChPT to give an acceptable description of our data. All the 
additional data introduced in this work satisfies this cut and we include all the data for pions with 
valence masses ra x , m y < 0.01 on the two 24 3 ensembles and all data for pions with valence masses 
mx, m y < 0.008 for the three 32 3 ensembles. For kaons we include all the valence light-quark 
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Parameter 


No FV Corrections 


With FV Corrections 


B 


4.12(7)GeV 


4.03(7) GeV 


f 


0.110(2) GeV 


0.112(2) GeV 


J 


0.05(7) GeV 2 


0.04(7) GeV 2 




-0.00000(7) 


-0.00005(7) 


j 


0.00050(5) 


0.00047(5) 


L (2) 

O 


-0.00003(4) 


-0.00005(4) 


L (2) 

8 


0.00055(2) 


0.00059(2) 




0.4856(4) GeV 


0.4854(4) GeV 


f(K) 


0.141(3) GeV 


0.143(3) GeV 


C f(K) 
J ' 


0.01(6) GeV 2 


0.01(6) GeV 2 


h 


0.0043(9) 


0.0046(10) 




0.023(1) 


0.024(1) 


h 


-0.0018(9) 


-0.0016(10) 




0.0058(2) 


0.0057(2) 




1.666(2) GeV 


1.666(2) GeV 




0.20(6) GeV" 2 


0.20(6) GeV" 2 



TABLE XXVII: Parameters of the global fit to our ensembles using NLO ChPT without finite-volume 
collections (second column) and with finite-volume corrections (third column). For the unitary theory the 
parameters are denned in Sect. IV Bl and for the partially quenched theory in appendix B of Ref. |l|]. 



masses in the above range for each fixed strange-quark mass. For this infinite- volume SU(2) NLO 
global fit the fitted parameters are presented in the second column of table lXXVIII The # 2 /dof 
for all the fits discussed here are given in table lXXVIII! We also perform the corresponding fits 
using the finite- volume chiral logarithm composed of a sum of Bessel functions [44]; resummed 
expressions are not available for our partially quenched fits. The parameters of the fit are presented 
in the third column of table lXXVIII In terms of the conventional LECs I3 and I4 the results are 



h = 2.82(16), 
h = 2.57(18), 



U = 3.76(9) 
U = 3.83(9) 



(Infinite Volume ChPT) 
(Finite Volume ChPT) . 



(59) 
(60) 
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Ansatz # 2 /dof 
NLO 0.72(46) 
NLO-fv 1.07(47) 
Analytic 0.60(44) 

TABLE XXVIII: Fit ansatze, mass ranges and uncorrelated # 2 /dof obtained in our analyses. The fits were 
performed for pion masses less than 420 MeV. 



Parameter 



Value 



c 



c 



-0.001(1) GeV 2 
7.45(9) GeV 
0.43(8) GeV 
0.123(2) GeV 
0.04(7) GeV 2 
0.85(2) 
0.56(9) 
0.2353(8) GeV 2 



Parameter 



Value 



Cf* 
c 2 



G 



C 



c 



Ik 
Jk 



C 



c 2 



3.67(4) GeV 
0.7(1) GeV 
0.149(2) GeV 
0.02(6) GeV 2 
0.34(1) 
0.52(10) 
1.666(2) GeV 
2.7(9) 



TABLE XXIX: Parameters of the global fit to our ensembles using the analytic ansatz. The parameters are 
defined in Eqs. (|49]> - d55l ). 



In table IXXIXI we present the parameters of the fit with the analytic ansatz over the same mass 
range as for the fits using SU(2) chiral perturbation theory, as explained in the previous paragraph. 
We find that analytic fits including a larger range of pseudoscalar masses give an acceptable un- 
correlated % 2 /&of but then the lightest data points were consistently missed by the fit by about 
one standard deviation. The utility of such extended fits for extrapolating to the physical point 
was therefore compromised and we therefore decided to restrict the range of masses used in the 
analytic fits. 

The global fit to many ensembles of partially quenched data is naturally a high dimensional space 
and so the exposition of the fits is best performed by looking at portions of the data in turn. In 
order to illustrate the quality of the fits, in the following subsections we display the fit and data for 
each physical quantity in turn. In total we have analysed five ensembles at two lattice spacings, 
and each ensemble has measurements at many partially quenched valence-quark masses. As it is 
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only feasible to present a subset of possible plots, in the following we display the dependence of 
each quantity on the valence quark masses at the lightest sea-quark mass (m/ = 0.005 for the 24 3 
ensembles and m/ = 0.004 on the 32 3 ensembles). The exception of course, is the mass of the 
Omega baryon mhhh which does not depend on the light valence-quark masses. We also display 
the unitary subset of data on both lattice spacings along with the mass dependence we infer from 
our fits in the unitary continuum limit. 

Before discussing the chiral and continuum behaviour of hadronic masses and decay constants in 
detail, we present in table lXXXl our results for the unrenormalised physical quark masses and the 
lattice spacings obtained from the three fits. In this table the quark masses are given in lattice 
units. The non-perturbative renormalization of the masses will be discussed in Sec .|VT1 where the 
values of the renormalized quark masses in the MS scheme will be presented. 





NLO 


NLO fv 


Analytic 


m/(32 3 ) 


0.00100(3) 


0.00102(3) 


0.00105(6) 


m,.(32 3 ) 


0.0280(7) 


0.0280(7) 


0.0279(7) 


fl- ! (32 3 ) 


2.280(28) GeV 


2.281(28) GeV 


2.282(28) GeV 


m,(24 3 ) 


0.00134(4) 


0.00136(4) 


0.00141(9) 


»\(24 3 ) 


0.0379(11) 


0.0379(11) 


0.0378(11) 


a- l (24 3 ) 


1.729(25) GeV 


1.729(25) GeV 


1.730(25) GeV 



TABLE XXX: Unrenormalised physical quark masses in lattice units and the values of the inverse lattice 
spacing a~ l for the 32 and 24 3 ensembles. 

1. Chiral and continuum behaviour of the Q.-baryon 

The Q. mass is fitted using Eq. (l45l) (or equivalently (l55l) ). The fit form for the Q. baryon does 
not change between the different ansatze and only very small differences arise from the different 
estimates of physical quark masses and hence of the lattice spacings. For illustration, Figure 1271 
shows the extrapolation of the Q. mass using the analytic ansatz. 
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1.78 
1.76 
> 174 

o 
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1.68 
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0.005 0.01 0.015 0.02 

nij (GeV) 

FIG. 27: The fit to the light-quark mass behaviour of the £2-baryon in the continuum limit obtained using 
the analytic ansatz. The corresponding plots using the infinite and finite-volume SU(2) ChPT ansatz are 
almost indistinguishable, differing only slightly in the estimates of the physical quark masses and the lattice 
spacings. 

2. Chiral and continuum behaviour of the pion mass 

We display the fits of the partially quenched pion masses using infinite volume NLO SU(2) par- 
tially quenched ChPT (i.e. to the partially quenched generalization of Eq. (1381) given in Eq. (B.32) 
of ref. lit) in figure 128] for the lightest 24 3 and 32 3 ensembles. As discussed in section IVCl we 
divide by the average valence-quark mass with the intention of enhancing the visibility of chiral 
logarithms. Figure [29] displays the corresponding fit of the same data but including finite-volume 
corrections. 

It is apparent that the infinite volume and finite volume NLO fits diverge rapidly from our data at 
larger masses, and this indeed is the reason why we were compelled to introduce the upper cut-off 
of 420 MeV for this analysis Jl]. 

We now consider the chiral extrapolation of the pion mass using the analytic form of Eq. (|49l which 
is shown in Fig. [30] Comparing Figs. 1281 and |29l with Fig. [30] suggests that data at substantially 
larger masses can be described by the analytic expansion, without any curvature terms in the 
ansatz. The division by the average valence quark mass in the plots, coupled to allowing the 
tangent not to pass through the origin (i.e. that the extrapolated rn\ at m x = m y = may not be 
equal to zero) allows the analytic fit to reproduce a structure that might otherwise be attributed to 
chiral logarithms. 

We emphasize that admitting the possibility that the constant term Cq k ^ allows for a pole in 
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FIG. 28: Global fits obtained using infinite volume NLO SU(2) chiral perturbation theory for the pion mass. 
The top-left panel includes the partially quenched data from the mi = 0.005 ensemble on the 24 3 lattice and 
the data points in the top-right panel are from the mi = 0.004 ensemble from the 32 3 lattice. In each case 
the curves correspond to the appropriate value of the lattice spacing. The points marked by the circles were 
included in the fit, whereas those marked by the diamonds were not. In the bottom two panels we zoom into 
the low-mass region, illustrating the fits to the points which were included (24 3 points on the left and 32 3 
points on the right). (For fixed m x , m y decreases as (am 9 ) 2 /m avg increases.) 



figure |30] in the unitary chiral limit. In fact we find that Cq* is numerically small and consistent 
with zero, C™ % = — 0.001 (l)GeV 2 . We stress again that while Goldstone's theorem implies the 
vanishing of the pion mass in the SU(2) chiral limit, this does not necessarily imply that Cq k = 0. 
Our model is that the linear ansatz is valid in the region between that where we have data and 
the physical point, and that if ^ then it is the curvature due to chiral logarithms below the 
physical pion mass which will force the pion mass to zero in the chiral limit. Nevertheless, from 
the fits we found that Cq* is consistent with zero. This is illustrated by the flat behaviour (within 
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FIG. 29: Global fits for the pion mass obtained using NLO SU(2) chiral perturbation theory with finite- 
volume corrections. In this case we only include the points which were included in the fit (mi = 0.005, 
24 3 points on the left and mi = 0.004, 32 3 points on the right) since the finite-volume corrections at larger 
masses are small. (For fixed m x , m y decreases as (am^) 2 /m avg increases.) 

the statistical precision) for the chiral behaviour of the unitary points for m^/mi in the continuum 
limit shown in the right panel in Fig. [31] Allowing for a non-zero value of Cq k does however 
lead to an amplified error for m|/m/ at the physical point. The left panel of Fig. [31] shows the 
corresponding plots for the infinite and finite-volume ChPT fits. 

Goldstone's theorem equally applies at vanishing valence-quark mass (m x = m y = 0) but with a 
non-zero sea-quark mass (m/ > 0). In contrast with the unitary case discussed in the previous 
paragraph where Cq % was consistent with zero, in the partially quenched direction we find that the 
corresponding constant Cq" -\-C^ % mi is non-zero, specifically C™ 71 = 0.43(8)GeV. This value for 
C™ % is much larger than might be created by propagating the mass dependence in m ; ies (m) through 
the term involving C" !?t ; the greatest mass dependence in m' K& occurs on our 24 3 ensembles in the 
partially quenched direction, but can at most generate a 1% correction to m and produces a term 
much smaller than the measured C^ 1 *. Further, the residual chiral symmetry breaking is four times 
smaller for the 32 3 ensemble which is also included in the global fit. Our results from this global 
analytic fit therefore require a curvature, most likely from partially-quenched chiral logarithms 
which are known to be larger than in the unitary direction, in order for Goldstone's theorem to be 
satisfied. 

It is also worth emphasizing that the discovery of chiral logarithms in lattice data from plots such 
as those in Figs.[28]to[30]is to a certain extent artificial. Inconsistency with LO chiral perturbation 
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FIG. 30: Global fit curves obtained using the analytic fit ansatz (|49l ) overlaying the simulated pion masses 
on the mi = 0.005, 24 3 ensemble (top-left) and the m/ = 0.004, 32 3 ensemble (top-right). Points marked by 
circles were included in the fit, those marked by diamonds were not. The simple linear expansion replicates 
the entire range of lattice data reasonably well with the description being rather better than NLO chiral 
perturbation theory at our larger masses. In the bottom two panels we zoom into the low-mass region, 
illustrating the fits to the points which were included (24 3 points on the left and 32 3 points on the right). 
(For fixed m x , m y decreases as (am^) /m avg increases.) 

theory is certainly indicated. Our linear fits suggest that the transformations made in displaying 
the data render even conclusions of genuine curvature, let alone unambiguous demonstration of 
logarithmic mass dependence, to be somewhat optimistic. In order to prove logarithmic behaviour, 
one should really change quark masses substantially on a logarithmic scale; our present lattice data 
supports only the weaker claim of consistency with logarithmic behaviour in the partially quenched 
direction. 
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m, (GeV) m, (GeV) 

FIG. 3 1 : Left panel: Pion mass fit for the SU(2) NLO fit form in the continuum limit, both with and without 
finite volume logarithms. We adjust the data points to the continuum limit using the a 2 dependence in our 
fit form and overlay these. Right panel: Chiral extrapolation of the pion mass using the analytic (l52l and 
infinite-volume NLO ChPT ansatze. 

3. Chiral and continuum behaviour of the pion decay constant 

We now turn to the chiral behaviour of f % and the extrapolation to the physical point. The leading 
term in all the fits contains an a 2 correction and we display the fits performed at non-zero lattice 
spacing combined with the unmodified lattice data and also our continuum predictions combined 
with the lattice data extrapolated to the continuum limit using the results of the fits. 
We display our fits obtained using infinite volume NLO SU(2) partially-quenched ChPT in Fig- 
ure[32l The corresponding fits including finite- volume corrections are shown in Figure[33l Finally 
Figure[34]displays the fits obtained using our analytic ansatz. Having performed the fits, we adjust 
our unitary data to the continuum limit using the fitting functions with the determined parameters 
and display the adjusted data in Fig. [35] together with the finite and infinite-volume NLO SU(2) 
ChPT fits (left panel) and the analytic fit (right panel). The effect of the adjustment to the con- 
tinuum limit is illustrated in Figurel36l where the fits are superimposed on the unadjusted unitary 
data. It can be seen from Figs.[35]and[36]that the adjustment to the continuum limit for the pion 
decay constant is very small. 

The predictions for f n extrapolated to the physical quark masses for each of the fits is given in 
table IXXXIl We anticipate the discussion of the global fits for fx which are presented in Sec IV E 61 
and mention that the predictions for fx extrapolated to the physical quark masses are given in 
table IXXXIIl and the predictions for fx/fn extrapolated to the physical quark masses are given in 
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m m 

FIG. 32: Global fits to the lattice data for the pion decay constant obtained using infinite- volume NLO 
SU(2) chiral perturbation theory. The top-left and top-right panels correspond to the 24 3 , m; = 0.005 and 
32 3 , mi = 0.004 ensembles respectively. Points marked by circles are included in the fits, while those with 
heavier masses marked by diamonds are not. In the bottom two panels we zoom into the low-mass region, 
illustrating the fits to the points which were included (24 3 points on the left and 32 3 points on the right). 
(For fixed m x , m y increases as afxy increases.) 

table muni 

We find that the NLO SU(2) fits underestimate the physical value at our simulated lattice spacings, 
and that this discrepancy is amplified a little by the extrapolation to the continuum limit. At each 
of our two lattice spacings, the analytic ansatz extrapolates close to the physical value of f n , but, 
with our ansatz for the form of the a 2 effects, the result becomes statistically inconsistent in the 
continuum limit. 

From the above discussion we see that using NLO ChPT to perform the chiral extrapolation for 
results in a value which is significantly smaller than the physical one. We recall that only data 
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0.012 0.014 




FIG. 33: Global fits to the lattice data for the pion decay constant obtained using NLO SU(2) chiral 
perturbation theory with finite-volume corrections. In this case we only include the points which were 
included in the fit (mi = 0.005, 24 3 points on the left and m/ = 0.004, 32 3 points on the right) since the 
finite- volume corrections at larger masses are small. (For fixed m x , m y increases as afxy increases.) 





NLO 


NLO fv 


Analytic 


J n 


0.121(2) 


0.123(2) 


0.128(2) 


J n 


0.120(2) 


0.122(2) 


0.127(2) 


/■continuum 
J 7t 


0.119(2) 


0.121(2) 


0.126(2) 



TABLE XXXI: Predictions for f n in GeV for each global fit ansatz at each simulated lattice spacing and in 
the continuum limit. 

limited to m K < 420 MeV was used in the analysis and note that the fits were performed using the 
chiral expansion with /, the decay constant in the SU(2) chiral limit, included in the expansion 
parameter Xi/i^xf) 2 - The downward curvature at low masses seen in Figure[35]can, of course, be 
reduced by replacing the mass-independent / by an artificial larger parameter such as the physical 





NLO 


NLO fv 


Analytic 


/■24 3 
JK 

f 32 3 
JK 

/■continuum 
JK 


0.147(2) 
0.147(2) 
0.146(2) 


0.148(2) 
0.148(2) 
0.147(2) 


0.152(2) 
0.151(2) 
0.151(2) 



TABLE XXXII: Predictions for in GeV for each global fit ansatz at each simulated lattice spacing and 
in the continuum limit. 
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FIG. 34: Global fits to the lattice data for the pion decay constant obtained using the analytic ansatz in 
Eq. (l52l) . The top-left and top-right panels correspond to the 24 3 , mi = 0.005 and 32 3 , m/ = 0.004 ensembles 
respectively. Points marked by circles are included in the fits, while those with heavier masses marked by 
diamonds are not. In the bottom two panels we zoom into the low-mass region, illustrating the fits to the 
points which were included (24 3 points on the left and 32 3 points on the right). (For fixed m x , m y increases 
as afxy increases.) 



fn or fu(fhi) measured at each quark mass used in the simulation. The curvature can also be 
partially absorbed by using a subset of terms that arise at NNLO. We have experimented with 
NNLO fits Ii46ll but find that the low-energy constants are insufficiently constrained by our data to 
be of practical use. Thus the resulting predictions for the physical value of f n depend strongly on 
the model assumptions used at NNLO. 

The observed O(10%) deviation found using NLO chiral perturbation theory is broadly consistent 
with the size of NNLO terms one might expect to be present at masses in the region of our data. 
Our data for f K vary from about 20% to 40% above the value of / obtained from our extrapolations 
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FIG. 35: Unitary data for f n adjusted to the continuum limit using each of the fit ansatze. The left panel 
compares the infinite volume and finite volume forms of the NLO SU(2) fit, while the right panel com- 
pares the analytic fit to the infinite volume NLO SU(2) fit. The horizontal solid line indicates the value 
/ w -= 130.4 MeV (the authors of ref. Q45D quote f K - = ( 130.4 ±0.04 ±0.2) Me V). 





NLO 


NLO fv 


Analytic 


{hlfn) li? 


1.216(9) 


1.205(9) 


1.184(9) 


(h/fn? 2 ' 


1.221(6) 


1.209(6) 


1.188(6) 


If If \ continuum 
{JK/Jn) 


1.229(8) 


1.215(7) 


1.194(7) 



TABLE XXXIII: Predictions for fx/ fx for each global fit ansatz at each simulated lattice spacing and in 
the continuum limit. 

and the square of these terms can be taken as being indicative of the expected NNLO terms. We 
might therefore expect them to be around 5-15% within our simulated mass range. 
The discrepancy of the prediction for the physical value of f K from the analytic fits is smaller than 
that found with NLO ChPT, but is nevertheless visible. The results at each of the two lattice spac- 
ings are statistically consistent with f % but lead to an underestimate in the continuum limit. Given 
the sign of the chiral logarithms at NLO, one might expect a linear ansatz to over-estimate rather 
than underestimate the prediction for the physical value. It is nevertheless striking that one cannot 
admit any significant non-linearity in this extrapolation and retain consistency with the physical 
value for f K . The simple analytic form used here appears to be a successful phenomenological 
model which is simpler and has fewer parameters than approaches based on ChPT with arbitrarily 
chosen analytic subsets of NNLO and NNNLO terms. 



76 



0.16 - 



0.15 



0.14 



> 

a 



! 0.13 



0.12 



0.1 1 




0.005 



LO analytic 
NLO SU(2) ChPT 
LO analytic 32" 
LO analytic 24 3 
NLO SU(2) ChPT 32 J 
NLO SU(2) ChPT 24 3 
32 3 data 
24 3 data 



0.01 0.015 
m. (GeV) 



0.02 



FIG. 36: Chiral extrapolation of the pion decay constant using the analytic d52l ) and ChPT (1421 fit ansatze. 
Here, the lattice results from the 24 3 and 32 3 ensembles are shown along with the mass dependence we 
infer both at each lattice spacing and in the continuum limit. The consistency of the two ensembles with 
each other and with this continuum limit is indicative of the size of lattice artefacts. The horizontal solid 



line indicates the value f n - = (130.4 ± 0.04 ± 0.2) MeV |4J]. 



It is of interest to pose the scientific question whether any of the fit ansatze could in principal be 
consistent with the experimentally measured pion decay constant? To answer this question we 
update the analysis of Ref. Q47D and include an artificially created data point for each ensemble 
that represents the experimental result in the continuum limit but includes our fitted a 1 correction 
at each non-zero lattice spacing. This is displayed in figure [37] and we find that the analytic 
ansatze could be consistent with an uncorrected ^ 2 /dof = 1.9(7), while NLO ChPT would fail 
to simultaneously fit our data and the physical point, with % 2 /&of =6(1) (infinite volume) and 
# 2 /dof =5(1) (finite volume). 

Of course, improved statistical errors, simulations at a third lattice spacing and larger physical 
volumes would give us better control of the continuum extrapolation and finite-volume effects. 
However, our main conclusion is that it is imperative to simulate with masses substantially nearer 
to the physical point; this will constrain both fit forms to give more consistent predictions. Ul- 
timately simulations will be performed directly at physical quark masses and will eliminate this 
error completely. We are currently generating new ensembles with a coarser lattice spacing, with a 
substantially larger volume and with very much lighter pion masses (for a preliminary discussion 
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FIG. 37: An artificial data point (the left-most data point in each panel) corresponding to the physical value 

n 

°f fx H45TI . but including our uncertainties in the lattice spacing, is added to the data for the pion decay 
constant from the five ensembles. The left-hand panel corresponds to the NLO SU(2) ChPT fits and the 
right-hand panel to the analytic ansatz. 



of these configurations see Ref. H&H ) precisely to address this issue. 

As an estimate of the systematic uncertainties in physical quantities we take the difference be- 
tween the results obtained using linear and finite- volume NLO ChPT analyses. This allows for the 
possible validity of the full NLO non-analyticity in the region of masses between the data and the 
physical point but also recognises that part of this extrapolation may be outside the range of valid- 
ity of NLO ChPT as suggested by the observation that the present data is surprisingly consistent 
with linear behaviour. Guided by the results for f K discussed above, we take as our central values 
for phenomenological predictions the average of the results obtained from our finite- volume NLO 
ChPT fits and our analytic fits. 



4. Chiral and continuum behaviour of the mass of the kaon 

We display our fits using infinite volume NLO SU(2) partially quenched ChPT in figure[38l Fig- 
ure[39]displays the corresponding fits of the same data with the finite- volume corrections included, 
while the analytic fits are displayed in figure [401 The corresponding unitary view of the data in the 
continuum limit is shown in figure EU All these plots are for results at the physical sea strange 
quark mass. 
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FIG. 38: Dependence of the kaon mass on the mass of the light valence quark with fits performed using 
infinite- volume NLO partially-quenched ChPT. The left panel shows the results from the 24 3 , mi = 0.005 
ensemble and the right panel from the 32 3 , mi = 0.004 ensemble. In each case the results are for the physical 
strange-quark mass. 
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FIG. 39: Dependence of the kaon mass on the mass of the light valence quark with fits performed using 
finite- volume NLO partially-quenched ChPT. The left panel shows the results from the 24 3 , mi = 0.005 
ensemble and the right panel from the 32 3 , mi = 0.004 ensemble. In each case the results are for the 
physical strange-quark mass. 

5. Chiral and continuum behaviour of /k 

We next discuss f%, the decay constant of the kaon. We display our fits using infinite-volume 
NLO SU(2) partially quenched ChPT in Figure|42l The following two figures display fits of the 
same partially quenched data to ChPT with finite-volume corrections (Figure|43]) and to the global 
analytic fit ansatz (Figure[44]). The NLO ChPT fit ansatze, both with and without finite-volume 
logarithms, are displayed for the unitary data adjusted to the continuum limit in figure |45| 
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FIG. 40: Dependence of the kaon mass on the mass of the light valence quark with fits performed using the 
analytic fit ansatz. The left panel shows the results from the 24 3 , m/ = 0.005 ensemble and the right panel 
from the 32 3 , m; = 0.004 ensemble. In each case the results are for the physical strange quark mass. 
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FIG. 41: Chiral extrapolation of the kaon mass using unitary data points adjusted to the continuum limit by 
the fitting ansatze. Here we compare results obtained using the infinite-volume NLO ChPT ansatz to that 
using finite volume logarithms (left panel) and to the analytic ansatz (right panel). 

The two panels in Figure |46] display the chiral behaviour of the actual unitary data from the two 
sets of ensembles (left panel) as well as of the data adjusted to the continuum limit (right panel). 
From these fits our final predictions for f% are given in table IXXXIIl and the corresponding results 
for 4 in table IXXXIIIl 



6. Predictions 



We now present our results for f % , fx and their ratio as well as for the physical bare quark masses. 
As discussed above, our central value for any physical quantity is taken to be the average of the 
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FIG. 42: Dependence of the kaon decay constant on the mass of the light valence quark with fits performed 
using infinite- volume partially quenched NLO ChPT. The left panel shows the results from the 24 3 , mi = 
0.005 ensemble and the right panel from the 32 3 , mi = 0.004 ensemble. In each case the results are for the 
physical strange quark mass. 



results obtained from analyses using the NLO SU(2) ChPT fit with finite volume corrections and 
those from the analytic fit. The difference between the analytic and finite-volume NLO SU(2) fits 
is taken as a systematic error. This procedure includes a NLO finite-volume correction, estimated 
from the difference between results obtained using NLO ChPT at infinite and finite volumes, and 
which is much smaller than the total systematic error here. 

Our predictions for pseudoscalar decay constants therefore contain systematic errors for finite 
volume effects, the chiral extrapolation, and residual chiral symmetry breaking, while the discreti- 
sation error is included indirectly by the fitting procedure: 



124(2) (5) MeV (61) 
149(2) (4) MeV (62) 
1.204(7)(25), (63) 

where we display the statistical and systematic errors separately. We note that the known, exper- 
imental value of f K influenced our choice to take the central value of physical quantities as the 
average of the results from the analytic and finite-volume NLO ChPT ansatze. The prediction for 
f K cannot therefore be considered unbiased, however as our aim is to select the most likely central 
value for phenomenologically important quantities such as fx/ fit and Bk our procedure is both 
appropriate and contains a prudent systematic error. 
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FIG. 43: Dependence of the kaon decay constant on the mass of the light valence quark. The left panel 
shows the results from the 24 3 , mi = 0.005 ensemble and the right panel from the 32 3 , m/ = 0.004 ensemble. 
In each case the results are for the physical strange quark mass. There are two curves plotted. The orange 
curve is the result one infers for the infinite volume, while the red curve is the result we obtain on the finite 
volume. As we do not adjust our data for finite volume effects, the red curve should go through our data. 
The orange curve also goes through our data which is an indication that the finite volume effects in our 
data are substatistical, and the difference between the orange and red curves at lighter masses indicates that 
one should expect substantial finite volume effects if one were to simulate at these lighter masses without 
changing our present volume. 

Applying the same procedure to obtain predictions for the physical bare quark masses for the 
/3 = 2.25 32 3 ensembles, we find: 

m Mrf = 2.35(8)(9)MeV and m s = 63.7(9) (l)MeV, (64) 

and these will be renormalised in the following section. The corresponding bare masses for the 
/3 = 2.13 24 3 ensembles can be obtained by dividing the results in (|64|) by the values of Z/ and Zf, 
in Table|XXVU 

7. Chiral and continuum behaviour ofr$ and r\ 

Finally in this section we apply the combined chiral/continuum extrapolation procedure to the 
scales tq and r\. Assuming a linear dependence for the light sea-quark mass dependence, and 
including a leading order a 1 term as before, the scales are independently fit to the form 

ri = c n + c rha a 2 + c rim ihi , (65) 



82 




FIG. 44: Dependence of the kaon decay constant on the mass of the light valence quark with fits performed 
using the analytic fit ansatz. The left panel shows the results from the 24 3 , mi = 0.005 ensemble and the 
right panel from the 32 3 , mi = 0.004 ensemble. In each case the results are for the physical strange quark 
mass. 
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FIG. 45: Chiral extrapolation of the kaon decay constant for unitary data in the continuum limit. We 
compare the NLO ChPT ansatz to the corresponding ansatz with finite-volume logarithms. 

where i = 0,1. Prior to the fit, the data are linearly interpolated to each of the physical strange 
quark masses obtained from the global fits and presented in Table lXXXl and the fit and the subse- 
quent extrapolation are performed using the corresponding physical light-quark mass and lattice 
spacings. 

The parameters and ^ 2 /d.o.f of the fits are given in Tables IXXXIVI and IXXXVI respectively, and 
plots showing the fits overlaying the data in the continuum limit are shown in figure |47J The fits 
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FIG. 46: Chiral extrapolation of the kaon decay constant for unitary data in the continuum limit. We 
compare the NLO ChPT ansatz to the analytic ansatz. The left panel displays the data and fits at non-zero 
lattice spacing, while the right panel displays the predicted results and correspondingly adjusted data points 
for the continuum limit. 

(a) r 



Parameter 


ChPT 


ChPT-fv 


Analytic 


c ro 


2.468(41) GeV- 1 


2.468(41) GeV" 1 


2.467(41) GeV" 1 


c r ,a 


-0.25(14) GeV 


-0.25(14) GeV 


-0.25(14) GeV 


c r ,mi 


0.42(1.23) GeV" 2 


0.44(1.23) GeV- 2 


0.47(1.23) GeV" 2 






(b) n 




Parameter 


ChPT 


ChPT-fv 


Analytic 




1.694(29) GeV- 1 


1.694(29) GeV" 1 


1.693(29) GeV" 1 


Cr\ ,a 


-0.15(11) GeV 


-0.15(11) GeV 


-0.15(12) GeV 


Cr\ ,rri[ 


-1.76(64) GeV" 2 


-1.76(64) GeV" 2 


-1.76(64) GeV" 2 



TABLE XXXIV: Parameters of the chiral/continuum fits to ro and r\. 

to ro appear to describe the data well by eye, and have a reasonable (uncorrelated) ^ 2 /d.o.f for 
the central value, but with a large deviation across the superjackknife distribution. The fits to r\ 
also appear to describe the data reasonably well, although there does seem to be a tension with the 
heaviest point on the 24 3 ensembles, which is likely responsible for the larger ^ 2 /d.o.f. As there 
are only five data points it is difficult to reach any stronger conclusions regarding the data: more 
ensembles and better statistics are needed. For the purpose of quoting a final result, we apply a 



PDG scale factor of ^/^ 2 /d.o.f to the statistical errors on each of the results. In order to retain 
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Quantity 



ChPT ChPT-fv Analytic 



1.35(1.66) 1.34(1.65) 1.31(1.63) 
2.69(2.39) 2.68(2.38) 2.66(2.37) 



TABLE XXXV: £ 2 /d.o.f of the chiral/continuum fits to r and r\. 



Quantity 


ChPT 


ChPT-fv 


Analytic 




2.469(39) GeV 


1 2.469(39) GeV 


1 2.468(39) GeV 1 


n 


1.690(29) GeV 


1 1.690(29) GeV" 


1 1.689(29) GeV 1 


n/ro 


0.6844(96) 


0.6844(97) 


0.6843(97) 



TABLE XXXVI: Continuum values of ro and r\ and the ratio n/ro at physical quark masses determined 
from a chiral/continuum fit using the lattice spacings and quark masses obtained from the global fits. 

the correlations between these quantities when the ratio is taken, the scale factor is applied to the 
difference of each jackknife sample from the mean. 

The continuum results for ro, r\ and their ratio at physical quark masses are given in table IXXXVII 
Using the procedure for combining the results obtained using the different chiral ansatze outlined 
in Section lVE3l and applying the PDG scale factor as above, gives: 

r = 2.468(45) stat (l) FV (l) z GeV- 1 = 0.4870(89) stat (2) FV (2) z fm, 

n = 1.689(47) stat (0) F v(l)^GeV- 1 = 0.3333(93) sta t(l)Fv(2)^fm, and (66) 

n/ro = 0.684(15) sta t(0)Fv(0) z , 

where the finite volume error arising from the different determinations of the lattice spacings and 
quark masses is smaller than the quoted precision on the ratio. % labels the error due to the chiral 
extrapolation. For comparison, the MILC collaboration recently obtained r\ = 0.31 17(6) ( jl^i) f m 
(~ 1.580(3)(^ 6 )GeV- 1 ) 049J], and also n = 0.317(7)(3) fm J~ 1 .61 (4) (2) GeV -1 ) and r = 
0.462(H)(4) fm (~ 2.34(6) (2) GeV -1 ) from an earlier study Q50Q . At this time we do not have 
an explanation of the discrepancy between our results in (l66l) and those of the MILC collabora- 
tion beyond noting the very different approaches to setting the scale and performing the chiral 
extrapolation. 
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FIG. 47: The scales ro (left) and r\ (right) corrected to the continuum limit, overlaid by the chiral/continuum 
fit. The extrapolated point at the physical light quark mass is shown as the grey cross. Here the lattice 
spacings and physical light quark mass were obtained from the global fits using the analytic ansatz. The 
fits using the quantities obtained with the ChPT and ChPT-fv global fit ansatze are almost indistinguishable 
from those shown in these figures. 

VI. LIGHT-QUARK MASSES 



The quark masses quoted in Eq. (l64~l) are the bare masses for the lattice action which we are using 
on the 32 3 ensembles with /3 = 2.25 corresponding to a lattice spacing a -1 ~ 2. 28 GeV. In order 
to be useful in phenomenological applications these results must be translated into renormalized 
masses in some standard continuum scheme. Therefore in Subsection lVI Al we determine the 
renormalization constants relating the bare masses in (|64l) to those renormalized in the MS scheme 
at a renormalization scale of 2 GeV In Subsection lVI Bl we then combine these renormalization 
constants with the bare masses in (l64~l) to obtain the renormalized masses, the LO LEC B MS (2 GeV) 
and the chiral condensate. 



A. Non-perturbative renormalization for quark masses 



The quark-mass renormalization factor which relates the lattice bare quark mass to that in the MS 
scheme is determined using non-perturbative renormalization (NPR) with the RI/SMOM schemes 
proposed in Ref. || 14f| as intermediate schemes. This is an extension of the Rome-Southampton 
NPR program in which the RI/MOM scheme was defined 115 ill . Quark masses renormalized in the 
RI/SMOM or RI/MOM schemes are obtained entirely non-perturbatively. Since it is not possible to 
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simulate in a non-integer number of dimensions, continuum perturbation theory is needed to match 
the results in either the RI/SMOM or the RI/MOM scheme and the target MS scheme. We stress 
however, that we completely avoid the use of lattice perturbation theory which often converges 
more slowly than continuum perturbation theory (PT). Since RI/MOM and any of the schemes 
proposed in jl4 \ are legitimate renormalization schemes, we exploit the freedom to choose an 
intermediate scheme to reduce its effect on the final result for the renormalized quark mass in the 
MS scheme and to have a better understanding of this uncertainty. 

Our earlier study fl 1 3M - used to normalize the quark mass on the 24 3 ensembles, applied the 
RI/MOM scheme to renormalize the quark masses and suffered from sizable systematic errors 
with two dominant sources. One of these is the truncation error in the perturbative continuum 
matching between the RI/MOM and MS schemes. This was estimated to be 6% for ji = 2GeV 
from the relative size of the highest-order term used (3 loop). The other is a non-perturbative effect 
arising because the strange quark mass is fixed close to its physical value, and the chiral limit is not 
taken for this quark. We estimated the corresponding systematic error on the quark-mass renormal- 
ization factor for a -1 = 1.73 GeV and /i = 2 GeV to be about 7%. As the strange-quark mass and 
the typical scale of spontaneous chiral symmetry breaking are almost the same, this error can be 
viewed as a general error due to contamination of non-perturbative effects (NPE). It was shown in 
Ref. [13] that changing the kinematics of momenta used to define the NPR scheme greatly reduces 
the contamination from unwanted non-perturbative effects and this will be discussed below. The 
actual implementation of the schemes with unconventional kinematics has been done in Ref. IU4I1 



carefully ensuring that the Ward-Takahashi chiral identities are satisfied. A pilot study Q52Q us- 
ing the new schemes demonstrated that it is a promising alternative to the conventional RI/MOM 
scheme with reduced systematic errors. In the present article we use two RI/SMOM schemes 
proposed in Ref. [14]. Preliminary results have been reviewed in Ref. [53]. 



An important technical improvement introduced since the previous study IU3I1 is the use of volume 



momentum sources for the quark propagators. This helps to reduce the statistical error greatly and 
in addition reduces the systematic error due to the dependence on the position of the local source 
used in []130. More details about the use of momentum sources can be found in Ref. [34j]. 
The mass renormalization factor Z m is conveniently calculated using the relation 



Z m = l/Zs= l/Zp, 



(67) 



where Z m , Z$, Zp are the quark mass, flavor non-singlet scalar and pseudoscalar renormalization 
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factors respectively. Here we are exploiting the important chiral symmetry properties of DWF. Our 
convention is that the renormalization factors multiply the bare quantities to yield renormalized 
ones: 

m R = Z m m 1 P a R =Z P P\ S a R =Z s S\ (68) 

where the left-hand sides are the renormalized mass, pseudoscalar and scalar densities and a is a 
flavour label, m in Equation (1681) is in physical units. The relations in Eq. (1671) are necessary for the 
Ward-Takahashi identities to hold for the renormalized operators. The RI/MOM renormalization 
condition on the amputated scalar vertex Els reads 

|*l T r[n 5 -/] = l. (69) 

Z q is the wave function renormalization factor, which can be determined using the trace condition 
on the local vector operator, 

|^[%r„] = i. (70) 

The vertex functions n depend on the incoming and outgoing momenta on the two fermion lines, 
T\(pin,p out ). The conventional RI/MOM scheme is defined using the forward vertex with pi n = 
Pout = P- The renormalization conditions Eqs. (l69l) , d70l) are applied by setting the renormalization 
scale /i to be the off-shell external momentum, \x 2 = p 2 , in the chiral limit. 
It is in principle possible to determine Z5 (= Zp) using the pseudoscalar vertex function instead of 
the scalar one in Eq. (l69l) . However, with the original RI/MOM choice for the external momenta, 
the pseudoscalar vertex couples to the zero-momentum pion, and the Green function diverges as 
l/m q as the quark mass m q — > at fixed p 



H54I1 . Therefore the pseudoscalar vertex cannot be 
used without some manipulation of the divergence (see e.g. 115511 ) and has not been considered in 
our previous publication II 1311 . This is in contrast with the RI/SMOM schemes described below 
which do not have such a pole as m q — > 0. Similarly, the axial-vector vertex can be used to 
determine Z q because Zy =Za- However, Z q obtained using the vector and axial- vector vertices at 
large but finite p 2 will differ because of the coupling of the axial current to the Goldstone boson 
115 ill . These differences are known to be of 0(1 / p 2 ) at high momentum from the operator product 
expansion 1 5 ll 15411 or from Weinberg's theorem of power counting for a Feynman diagram [13]. 
In Ref. II 1311 . the average of the vector and the axial- vector vertex was used to determine Z q and 
the difference was included in the systematic error, though the corresponding 1% error is sub- 
dominant. 
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The caveats mentioned in the two preceding paragraphs are both connected to the RI/MOM 
scheme and its channel with an "exceptional momentum"; specifically, the momentum transfer 
q = pi n — p out = 0. This is the reason for the large NPE error. It was demonstrated that the use 
of non-exceptional momenta pi„ — p out ^ reduces the NPE effect significantly. The RI/SMOM 
schemes are designed so that all channels have non-exceptional momenta. For quark bilinear 
operators we choose to have p\ n = p 2 out = q 2 and hence introduce the name "Symmetric Mom" 
(SMOM) schemes. The two schemes RI/SMOM and RI/SMOM yf( are defined with this kine- 
matical choice but differ in the T-projection operators which are used to define the wave function 
renormalization. For the vector (axial- vector) vertex function the projector g'q^/q 2 {Jsi^ii /q 2 ) is 
used in the RI/SMOM scheme and y M (y 5 y M ) as in Eq. CZQl) is used for RI/SMOM yfi . The standard 
/ (Ys) spinor projector is used for the scalar (pseudoscalar) vertex in both new schemes. 
The conversion factors from the RI/SMOM and RI/SMOM^ schemes to MS have been calculated 



at one-loop order in Ref . 111411 and recently to two-loop order 11151 . 



16] 



C m (RI/SMOM-»MS,/i) = 1- ^^0.646- (^^j (22.608 +4.014/1/) •■■ (71) 

C m (RI/SMOM^ — >-MS,jix) = 1- \^^jr) 1.979- f^j^J (55.032 + 6.162/1/) ••• (72) 

where the coefficients have been rounded to the third decimal place. Evaluating these factors at 
/i = 2 GeV we have 

C m (RI/SMOM^MS,jU = 2GeV,n/ = 3) = 1 -0.015 -0.006- •• , (73) 
C m (RI/SMOM r -> MS, ju = 2GeV, n f = 3) = 1 - 0.046 - 0.020 • • • . (74) 

In the RI/MOM and RI'/MOM schemes the conversion factors are known to three-loop order J^, 



51: 



C m (RI/MOM->-MS, J u = 2GeV,n/ = 3) = 1 -0.123 -0.070-0.048 + ••• , (75) 
C OT (Rl7MOM->MS, J u = 2GeV,n/ = 3) = 1-0.123-0.065-0.044 + ---. (76) 

We note that, at least up to two-loop order, the convergence of the series relating the new SMOM 
schemes to MS is considerably better than for the RI/MOM scheme. As already mentioned, the 
truncation error of the RI/MOM scheme was estimated from the size of the highest order term 
available (3 loop). Having in addition two intermediate SMOM schemes, we can expect to have a 
more reliable estimate of the truncation error. 
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We now turn to the numerical evaluation of the renormalization factors. At each value of /3, we 
use data obtained at the three light-quark masses: m; = 0.004, 0.006 and 0.008 for the finer 32 3 
lattice and m/ = 0.005, 0.01 and 0.02 for the coarser 24 3 lattice. 20 configurations were analyzed 
for each point. The ratio of quark wavefunction and local axial current renormalization factors is 
calculated from the average of vector and axial- vector vertex functions, 

|l = I(A v +Aa), (77) 
with projected and traced vertex functions: 

Ar MOM = _^_ Tr[n ^. | ^ ] ^ A™M = _^_ Tr[n ^. 75| ^ ]; (7g) 



for the RI/SMOM scheme. Here q^ in the continuum RI/SMOM scheme [|14l1 has been replaced 
with the q^ = sin^), as the derivative for the divergence of the current in the continuum theory 
is naturally replaced by the symmetric difference on the lattice. A remarkable feature of the 
RI/SMOM scheme is that in the chiral limit Ay = A^ holds non-perturbatively, in contrast to Ay ^ 
Aa for RI/MOM scheme due to spontaneous symmetry breaking (SSB). In principle there could 
still be a small difference for the lattice RI/SMOM scheme with non-zero m res , which, however, is 
negligible in the momentum range we use Q52Q . Using the continuum Ward-Takahashi identities, 
one can also show the equivalence of Z q in the RI/SMOM and RI'/MOM schemes Ill4n . 
The RI/SMOM^ scheme is defined using the conventional projectors, 

A RI/SMOM r „ 1 _ , , RI/SMOMy., 1 _ r „ , 

A v =^ n v»-y^ and A r '=^Tr[n AM . r5 r M ]. (79) 

Although these projectors are superficially the same as those used in the RI/MOM scheme, it 
should be remembered that the kinematics is different in the two cases with no exceptional chan- 
nels in the Green functions used to define the RI/SMOM^ scheme. 

The product of mass and wavefunction renormalization factors is calculated from the average of 
scalar and pseudoscalar vertex functions, 

1 

r 

with 



Z m Z q = -(A s + A P ), (80) 



A 5 = ^Tr[iVl] and A P = ^[U P ■ Js] , (81) 

again defined with the SMOM kinematics for the vertex functions. While A5 = Ap holds to all 
orders in perturbation theory with naive dimensional regularization, by using Weinberg's power- 
counting scheme we see that they can in general differ by terms of 0(1 / p 6 ) Jl3l1 . The difference 
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FIG. 48: Ap — As as a function of p 1 [GeV 2 ] for fine (32 3 ) and coarse (24 3 ) lattices. A straight line with 
\/p slope but arbitrary normalization is drawn to guide the eye. 

Ap — As after the chiral extrapolation is plotted in Fig. [48] as a function of p 2 (in physical units) 
for both the 24 3 and 32 3 lattices. The figure confirms the expected approximate 1/p 6 scaling. The 
unwanted non-perturbative effect from SSB is small and the introduction of non-exceptional mo- 
menta has had the expected effect. This is in contrast to the RI/MOM scheme with the exceptional 
channel, where the same difference behaves as 1/ (mp 2 ), and thus diverges in the chiral limit at 
finite p 2 . 

The mass renormalization factor Z^, with o = RI/SMOM or RI/SMOM^ , is given by combining 
Eqs. C77]) and ([80]), 

_ 1 A S + A P 

z '" " z7X*TXf (82) 

In calculating the ratio of vertex functions in Eq. (l82~l) we take the average of S and P or V and A 
for each light-quark mass and then fit with a quadratic (c+c'(m/ + m res ) 2 ) or linear c + c"(mi +w res ) 
formula to obtain the value c in the chiral limit for the numerator and denominator. For illustration, 
the extrapolation for the numerator using the quadratic formula is shown in Fig. |49l where the 
observed mass dependence is seen to be very small. Because of the very mild mass dependence, 
to the precision with which we quote our results and errors, the quadratic and linear extrapolation 
formulae lead to exactly the same quark-mass renormalization factor and error. Finally taking the 
ratio and combining with Zy gives the mass renormalization factor in the RI/SMOM schemes. The 
renormalization factor in the MS scheme at a scale /i = 2 GeV is obtained by first matching the 
scheme o to MS at /I 2 = pj n = p 2 ut = q 2 using Eqs. (1711 and (|72|) and then running to 2 GeV using 
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FIG. 49: Chiral extrapolation of (A/> + A^)/2 for the fine (32 3 ) lattice for each p 2 point. 
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FIG. 50: Z™° Mn ' (ju) and Z™ OM (jU) as functions of \i 2 = p 2 , and Z,^ s (2GeV) from the SMOM or 
SMOM 7fi schemes as function of matching scale squared p 2 for the fine lattice. The interpolation points are 
shown with the error bar at p 2 = (2 GeV) 2 . 



the three-loop anomalous dimension in the MS scheme. We use the four-loop QCD beta functions 
1 58] to calculate 0Cs 3 \n) for running and matching as shown in Appendix A of Ref. Jl3| . The 
relevant parameters taken from the 2008 Particle Data Group 04511 are 

ai 5 \m z ) = 0.1176, m z = 91.1876 GeV, m h = 4.20 GeV and m c = 1.27 GeV, (83) 

where the quark masses are in the MS scheme at the scale of the mass itself, e.g. = /«^ s (m^) . 



In Fig. [50] we plot Z, 



,SMOM 
m 



(H) and Z™ OM (ju) in the SU(2) chiral limit as functions of \l 2 = p 2 
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for the 32 3 ensembles. In addition we also plot Z,^ s (2GeV) as functions of the matching scale p 2 
obtained with SMOM and SMOMy M as the intermediate schemes. In an ideal situation, i.e. one in 
which the errors due to NPE contamination, truncation of perturbation theory and lattice artifacts 
are all small, the results obtained using the two intermediate schemes would give the same results 
for Z,^ s (2GeV), and the results would be independent of (pa) 2 . Since we have observed that the 
NPE error is small, the difference between the two sets of results is mostly due to the truncation of 
perturbation theory and lattice discretization errors. The observed decrease in this difference as p 2 
increases is consistent with the expected behaviour of the truncation error. Conversely, since the 
truncation error increases as p 2 decreases, taking the limit (pa) 2 — > 0, which is a typical treatment 
to eliminate the discretization error, is not an appropriate procedure. We therefore choose instead 
to evaluate Z m by taking an intermediate reference point p 2 = (2 GeV) 2 , for both the 24 3 and 32 3 
lattices. In this way, as we take the continuum limit of the renormalized quark mass, the leading 
(pa) 2 discretization error associated with the non-perturbative renormalization will be removed. 
There is a subtlety due to lattice artefacts which are not 0(4) invariant and which are responsible 
for the non-smooth (pa) 2 dependence in the figure. A term like a 2 Y t ^ i (Pn) 4 /p 2 , whose presence 
has been demonstrated in the conventional RI/MOM scheme for Wilson quarks [59], could exist 
also in the SMOM schemes. Such a term would manifest itself as scattered data around a smooth 
curve in p 2 , and the size of the scatter is expected to be comparable to the leading (pa) 2 error as 
both are of the same order in a 2 . This appears to be compatible to what is shown in the figure. 
Of course, it would be very helpful to know these terms, but in the absence of this knowledge we 



include this scatter in the systematic error by inflating the error by a factor A/z 2 /dof. The results 



are 



Z ; r (32) (M = 2GeV,n / = 3;sMOM r J = 1.573(2), (84) 
Z^ s(32) (Ai = 2GeV,n / = 3;sMOM) = 1.541(7). (85) 

The final arguments on the left-hand sides denote the choice of intermediate scheme. The error on 
the right-hand sides is the combination of the statistical fluctuations and the scatter of the points 
around the linear fit. The central values and errors are shown in the figure at the reference point, 
7? 2 = (2GeV) 2 . 

The 24 3 coarser lattice has been analyzed similarly for the mi = 0.005, 0.01 and 0.02 ensembles 
and the results are shown in Fig. [5TJ The mass renormalization factors on the 24 3 lattice for the 
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FIG. 51: Same figure as Fig. [501 but for the coarse 24 lattice. 



two intermediate SMOM schemes are: 

Z^ s(24) (/x = 2GeV,n / = 3;sMOM 7M ) = 1.578(2), (86) 
zjf (24) (/i =2GeV,n/ = 3;sMOM) = 1.534(10). (87) 

In Eq. (l64l) we have presented the bare quark masses for the fine 32 3 lattice and in Table IXXVII 
we give the ratios of equivalent bare masses on the 24 3 and 32 3 lattices. Because of the different 
0(a 2 ) artefacts for the light and heavy quark masses, there are two such ratios Z/ for the ud quarks 
and Z/ 2 for the s quark. These ratios Z/ and Z/ ? are also the scheme-independent ratios of the 
renormalization constants on the course and fine lattices. We now use these ratios to estimate 
the difference of the MS renormalized masses with the SMOM and SMOMy M schemes in the 
continuum limit. The continuum extrapolation of Z^ and Z^ /Zi or Zm^/Z/, will remove the 
(pa) 2 error in the non-perturbative renormalization. Thus, if a difference is found, it can largely be 
attributed to the truncation error of the perturbative matching. Performing such an extrapolation 
we find 

Z m/ S(32)C (M = 2GeV,n / = 3;sMOM yA1 ) = 1.527(6), (88) 
Z m/ S(32)C (M = 2GeV,n / = 3;sM OM ) = 1.511(22), (89) 

for the ud quark, and 

z 2 (32)C (i u=2GeV ^/ = 3 ; sMOM 7M) = 1-510(6), (90) 
Z ,^ (32)C (M = 2GeV,n/ = 3;sMOM) = 1.495(22) (91) 
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ensemble fine (32 2 ) coarse (24 3 ) coarse (16 3 ) ri31 

intermediate scheme RI/SMOM RI/SMOM RI/MOM 



PT truncation error 


2.1% 


2.1% 


6% 


m s / 


0.1% 


0.2% 


7% 


(A P -A 5 )/2 


0.5% 


0.6% 


N.A. (oo) 


(Aa-A v )/2 


0.0% 


0.0% 


1% 


total 


2.2% 


2.2% 


9% 



TABLE XXXVII: Systematic error budget for Z^ s (2GeV) with intermediate RI/SMOM schemes (this 



work) and RI/MOM scheme 111311 . 



for the s quark. Note that because these factors multiply m„^(32 3 )/a(32 3 ) or m 5 (32 3 )/a(32 3 ) 
presented in Eq. (|64|) to give the MS mass in the continuum limit, they are made to absorb the 
0(a 2 (32 3 )) discretization error in these bare quark masses on the fine lattice. Because of this, 
as well as the fact that the Z m 's are free from 0(a 2 ) errors originating from the SMOM non- 
perturbative renormalization, we have put additional suffix "c" as "continuum" to distinguish them 
from Zm S< ' 32 ' ) . The existence of a mass dependent contribution to the 0(a 2 ) artefacts gives rise to 
the different Z m for the light and heavy-quark masses. From the two different estimates of the MS 
renormalization factors with the SMOM and SMOMy^ intermediate non-perturbative schemes, 
we choose to take SMOMy M for our central value. The reason is that the scatter about the linear 
behaviour observed for the SMOM scheme in Figs.[50]and[5T]is much larger. Although the effect 
of the scatter has been taken into account in the error, we consider the continuum extrapolation 
from the SMOM scheme to be less reliable. The difference in the central values of ° in 



Eqs. (1881) and (1891) is about 1%, and this is also the case for the difference between the central 
values of Z^ 32 ^ in Eqs. (l90l) and (|91~I) . These differences of about 1% give an indication of the 
possible size of the truncation error of the perturbative two-loop matching to MS (it should be 
noted however, that the errors in the renormalization factors in the SMOM scheme are even a 
little larger). Another estimate of the truncation error of the matching is obtained by evaluating 
the size of the two-loop term in Eq. (1741) , resulting in 2.1% for the SMOM^ scheme. In order 
to be conservative, we shall take the latter as our estimate. Other systematic errors arise from the 
fact that the simulated strange mass is non-zero and from the small difference in the scalar and 
pseudoscalar vertices due to the residual spontaneous symmetry breaking effects. The first error 
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is estimated from the response of scalar and pseudoscalar vertex functions to the variation of the 



light-quark mass [13]. From the flat behaviour of Ap + A5 on the light-quark mass in Fig. [49] it can 



be seen that this uncertainty is small. The error estimates are compiled in Table IXXXVIII In the 



table, the corresponding errors from the RI/MOM analysis [13] are shown for comparison. All 
errors have become significantly smaller for the new SMOM schemes. Now our final values for 
the MS renormalization factor read 

z ,nf (32)C (M = 2GeV,n / = 3) = 1.527(6)(33), (92) 
Z 2P 2)C (M = 2 GeV ^/ = 3 ) = 1.510(6)(33), (93) 

where the first error is the statistical uncertainty inflated to take into account the scatter about the 
linear behaviour due to 0(A) non-invariant effects (as explained above) and the second is due to 
the remaining systematic effects and is dominated by the 2.1% truncation error of the perturba- 
tive matching. Here we have not taken into account the statistical fluctuation of Zy, which will be 
properly included in the calculation of the renormalized quark masses described in the next subsec- 
tion. The corresponding renormalization factor for the light-quark mass on the coarse 24 3 lattice 
is Z^ s(24)c (jU = 2 GeV,n/ = 3) = Z r Z™ {3Z)c (ii = 2 GeV,n f = 3) = 1.498 (6) (33). This value 
is consistent with our earlier estimate of the same quantity using RI/MOM as the intermediate 
scheme, 1.656(157) jl3\ \. but now with a considerably reduced error. 



B. Renormalized quark masses 

After the detailed discussion of the quark-mass renormalization, it is now straightforward to com- 
bine the renormalization constants in Eqs. (l92l) and (|931) with the physical bare quark masses on 
the 32 3 lattice in Eq. (|64l ) to obtain the light and strange quark masses renormalized in MS scheme: 

m^f(2GeV) = 2ff (32)c (/i = 2GeV,n / = 3) -m ud (32 3 ) •«- 1 (32 3 ) 

= 3.59(13) s t a t(14) sys (8) ren MeV, (94) 

mp(2GeV) = Z^ (32)c (/i = 2GeV,rc / = 3) -m,(32 3 ) - a - l (32 3 ) 

= 96.2(1.6) stat (0.2) sys (2.1) ren MeV, (95) 

where the three errors on the right-hand side correspond to the statistical uncertainty, the system- 
atic uncertainty due to the chiral extrapolation and finite volume, and the error in the renormaliza- 
tion factor. We recall that for the error due to the chiral extrapolation we conservatively take the 
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full difference of the results obtained using the finite-volume NLO SU(2) and analytic fits and for 
the central value we take the average of these results. We estimate the finite-volume effects from 
the difference of the results obtained using finite volume and infinite- volume NLO ChPT fits and 
combine these errors in quadrature. The finite-volume errors prove to be small. The error in the 
renormalization factor includes those in Eqs. (l92l) and (l93l) . 
The ratio of the s and ud quark masses is 

-^=26.8(0.8) stat (l.l) sys . (96) 
m ud 

We end this section by presenting our results for the leading-order LEC B and the chiral conden- 
sate. Using the finite- volume NLO ChPT fits we find 

5^(2GeV) = Z™ {32) -\ii = 2GeV,/ V = 3) -5(32 3 ) ■a- 1 (32 3 ) = 2.64(6) stat (6) sys (6) ren GeV. 

(97) 

Combining this result with the pion decay constant in the chiral limit, also obtained using the 
finite-volume NLO ChPT fits the chiral condensate is found to be 

[E M ^(2GeV)] 1 / 3 = [/ 2 J B(2GeV)/2] 1 / 3 = 256(5) stat (2) sys (2) ren MeV. (98) 

In Eqs. (|971) and ( |98T ) the second error is only due to finite volume corrections estimated from the 
difference of finite and infinite volume NLO ChPT fits. 



VII. TOPOLOGICAL SUSCEPTIBILITY 

The topological charge Q, defined on a single Euclidean space-time configuration, and its sus- 
ceptibility, Xq> are interesting quantities to calculate. While Q depends only indirectly on the 



quark masses, leading order SU(2) ChPT (60, 



611] predicts a strong dependence of Xq on the light 



sea quark mass with Xq vanishing linearly as m/ — > 0, suggesting that Xq ma y show important 

dynamical quark mass effects. 

In the continuum Q and Xq are defined by 

Q=^jd 4 xG^(x)G^ v (x) and X Q = (Q 2 )/V, (99) 

where V is the four- volume of the lattice, G^ v (x) is the gluon field strength tensor and G^ v {x), its 
dual. In the continuum, Q is integer valued and related to exact chiral zero modes of the massless 



Dirac operator by the Atiyah-Singer index theorem [62]. For sufficiently smooth gauge fields it 
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is possible to find a lattice expression which will always evaluate to an integer 116311 . as in the 
continuum limit. However, in the calculation reported here the necessary smoothness condition 
is not obeyed and we instead replace the right-hand side of Eq. (l99l) by a sum of Wilson loops 
chosen to approximate the G/ J , v (x)G lxv (x) product in Eq. (|99l) . Specifically we employ the "five- 
loop improved" (5Li) definition of the topological charge proposed in Ref. [64J] which at tree level 
is accurate through order a 4 . However, before evaluating this lattice expression for the topological 
charge, we smooth the links in the lattice by performing a series of APE smearing steps 1 6^1 



660 . The smearing parameter was set to 0.45, and 60 smearing sweeps were performed before 



measuring Q. The results are insensitive to the choice of these parameters. 

In Fig. [52] the Monte Carlo time history of Q is shown for each ensemble of gauge fields in our 
study. For each case, the update algorithm RHMC II [1] was used, except for the first 1455 
configurations for the m/ = 0.01 ensemble where the RHMC and RHMC I algorithms were 
used. In [|l|] it was shown that RHMC II is more effective in changing the gauge field topology, 
and therefore produces shorter auto-correlation times. The data for the first half (up to trajectory 
5000) of both 24 3 ensembles is repeated from fll]]. Figure [52] shows clearly the expected slowing 
of the rate of change of topological charge when moving towards the continuum Q6TD and, to a 
lesser degree, when decreasing the quark mass. The integrated auto-correlation times for Q for 
the smaller lattice spacing ensembles are shown in Fig. [2] While this figure is consistent with the 
autocorrelation times reaching a plateau of about 80 time units when integrated over an interval 
of about 200 time units, the exploding errors make this conclusion highly uncertain. Scanning 
Fig.[52]by eye, one might argue that the auto-correlations could be 500 time units, or longer. For 
example, note the large fluctuation to negative Q beginning around time unit 4750 for m/ = 0.006. 
The distributions of topological charge for each ensemble are shown in Fig. [53] The distributions 
become narrower as the quark mass is decreased. For the smaller lattice spacing, they also appear 
to exhibit non-Gaussian-like tails, or humps at large \Q\. 

Because of the parity symmetry of our calculation, the average of the pseudo-scalar quantity (Q) 
vanishes. However, %q remains non-zero and at leading order in SU(2) chiral perturbation the- 



ory i6 



6(1 



61Q is given by 



XQ = Z { ± + ±Y =E ^L, (10 0) 
\m u m d J m u + m d 

where E = Bf 2 /2 is the chiral condensate coming from a single flavor in the limit of vanishing up 
and down quark mass. 
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At one-loop in chiral perturbation theory Q68[] , 



1 1 ^ 1 



Xq = L — + — x 

(3 o JTl^ YYl III,/ \ 

[ -(w> mllos rf +K6(m " +mi)+2{2Kl+Ks) ^) t <101) 

= (l-^jmjlog^f + (2X 6 + 2X 7 + X 8 )m,), (102) 

where .K, = 128 EL,// 4 are proportional to the Gasser-Leutwyler NLO LEC's Q68Q , and in the last 
line the formula is evaluated for degenerate quarks. In contrast to other quantities considered in 
this paper, we do not attempt to characterize or evaluate the corrections to Eqs. (11011) or (11021) 
which come from non-zero lattice spacing. That interesting question is left for future work. 
In Tab. IXXXVIIII values of (Q) and Xq for each ensemble of configurations are summarized. To 
test for the expected auto-correlations, the data were blocked into bins of various sizes ranging 
from 10 to 600 time units. The quoted values of the statistical errors resulted when the block sizes 
were taken large enough that the errors no longer changed significantly. The block sizes are given 
in Tab. IXXXVIIII For all cases the first 1000 time units were discarded for thermalization. 
The dependence of Xq on the light quark mass is shown in Fig. [5H All of the data points lie 
above the LO curve (dashed line), all but the lightest significantly so. The result of the fit (^ 2 /dof 
Rj 13/4 fa 3) to the NLO formula Eq. (11021) is also shown. Since we have not determined K-j 
in Eq. (11021) from other means, we treat the linear combination of LEC's as a single, new, free 
parameter in the fit and find (2K^ + 2Kj + Kg) = 19.8(6.3). Except for the lightest data point, 
there is scant evidence for large 0(a 2 ) errors, though the statistical errors on the heavier two 
points with cT x = 2.284 are somewhat large. Omitting the former point in the fit leads to a more 
acceptable value of ^ 2 /dof w 1.5, suggesting the lightest point may be systematically low due to 
long auto-correlations in Q that are not well resolved in our finite Markov chain of configurations. 
Despite these limitations, the data appear to show a dependence on the light sea quark mass that is 
consistent with the dictates of NLO SU(2) ChPT 



VIII. CONCLUSIONS 

We have presented results from simulations using DWF and the Iwasaki gauge action for lattice 
QCD at two values of the lattice spacing (a _1 = 1.73 (3) GeV and a~ 1 =2.28 (3) GeV) and for uni- 
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TABLE XXXVIII: Topological charge and susceptibility. The measurement frequency, "meas. freq.", and 
"block size" are given in units of Monte Carlo time. 



mi 


meas. 


freq. block size 


(Q) (Q 2 ) X (GeV 4 ) 


0.005 


5 


50 


0.49 (25) 28.6 (1.4) 0.000290 (14) 


0.01 


5 


50 


-0.22 (37) 45.2 (2.5) 0.000458 (25) 


0.004 


4 


200 


0.59 (42) 11.4 (1.1) 0.000148 (14) 


0.006 


4 


200 


-0.07 (64) 24.8 (4.3) 0.000322 (55) 


0.008 


4 


400 


0.64 (100) 27.9 (5.6) 0.000363 (72) 
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FIG. 52: Monte Carlo time histories of the topological charge. The light sea quark mass increases from top 
to bottom, (0.005 and 0.01, 24 3 (top two panels), and 0.004-0.008, 32 3 ). Data for the 24 3 ensembles up to 
trajectory 5000 were reported originally in ||l|] and the results from the new ensembles are plotted in black. 
Most of the data was generated using the RHMC II algorithm (red and black lines). The RHMC (green 
line) and RHMC I (blue line) algorithms were used for trajectories up to 1455 for the mi = 0.01 ensemble. 
The small gap in the top panel represents missing measurements which are irrelevant since observables are 
always calculated starting from trajectory 1000. 
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FIG. 53: Topological charge distributions. Top: 32 3 , m, = 0.004-0.008, left to right. Bottom: 24 3 , 
mi =0.005 and 0.01. 

tary pion masses in the range 290-420 MeV (225-420 MeV for the partially quenched pions). The 
raw data obtained at each of the two values of /3 was presented in SectionslIIIl and ITVl respectively 
and the chiral behaviour of physical quantities on the 24 3 and 32 3 lattices separately was studied 
in Appendix|Aj The main aim of this paper however, was to combine the data obtained at the 
two values of the lattice spacing into global chiral-continuum fits in order to obtain results in the 
continuum limit and at physical quark masses and we explain our procedure in Section|V] In that 
section we define our scaling trajectory, explain how we match the parameters at the different 
lattice spacings so that they correspond to the same physics and discuss how we perform the ex- 
trapolations. We consider this discussion to be a significant component of this paper and believe 
that this will prove to be a good approach in future efforts to obtain physical results from lattice 
data. Although we apply the procedures to our data at two values of the lattice spacing, we stress 
that the discussion is more general and can be used with data from simulations at an arbitrary 
number of different values of /3. In the second half of Section|V]we then perform the combined 
continuum-chiral fits in order to obtain our physical results for the decay constants, physical bare 
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FIG. 54: Topological susceptibility (24 3 (squares), 32 3 (circles)). The dashed line is the prediction from LO 
SU(2) chiral perturbation theory (Eq. (1 1001 )) with the chiral condensate computed from the finite volume 
LEC's given in Table IXXVIII The solid line denotes the result of the single-parameter fit to the NLO 
formula given in Eq. (11021) . 

quark masses (which are renormalized in SectionlVIl) and for the quantities r$ and r\ defined from 
the heavy-quark potential. For the discussion below, it is important to recall that we use the phys- 
ical pion, kaon and Q. masses to determine the physical quark masses and the values of the lattice 
spacing and we then make predictions for other physical quantities. 

In contrast to most other current lattice methods, the DWF formulation gives our simulations 
good control over chiral symmetry, non-perturbative renormalization factors and flavor symmetry. 
This control allows us to measure and use, as either inputs or predictions: pseudoscalar decay 
constants, as well as their ratios; pseudoscalar masses; baryon masses; weak matrix elements 
and static potential values, limited only by the statistics achievable for these observables. The 
ability to predict many observables from the same simulations, provides evidence for the general 
reliability of the underlying methods. The good properties of DWF also allow us to test scaling, 
over this wide range of observables, at unphysical quark masses, since there are no flavor or chiral 
symmetry breaking effects to distort a test of scaling. We find scaling violations at the percent 
level, which supports including scaling corrections in only the leading order terms in our light- 
quark expansions. 
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As we reduce the quark masses used in the simulations, it is frustrating that there remains a doubt 
as to the best ansatz to use for the chiral extrapolation. We know of course that for sufficiently light 
u and d masses the behaviour is given by SU(2) ChPT; what we don't know is what "sufficiently 
light" means in practice. While in the range of quark masses accessible in our simulations, corre- 
sponding to 290 - 420 MeV for unitary pions and 225 - 420 MeV for partially quenched pions, our 
data are consistent with NLO SU(2) ChPT, we have seen that they are also consistent with a simple 
analytic ansatz leading to an inherent uncertainty in how best to perform the chiral extrapolation. 
This is particularly well illustrated in the study of f K , see Fig. [35] for example, where the data is 
well represented by all three ansatze (including NLO SU(2) ChPT with finite- volume corrections), 
but the extrapolated values differ as seen in Table|2QOO|/^ = 121(2) MeV from the NLO ChPT 
analysis with finite- volume corrections and f n =126(2) MeV using the analytic ansatz. Since a 
complete NNLO ChPT analysis is not possible with the available data, we have resisted the temp- 
tation to introduce model dependence by including only some of the higher order corrections and 
for our current "best" results we take the average of the two values and include the full difference 
in the systematic uncertainty obtaining f K = 124(2) (5) MeV. In Section lVE3l we investigated the 
increase in ^ 2 /dof if the fits are required to pass through the physical value 130.7(4) MeV up to 
corrections from lattice artefacts and found # 2 =1.9(7) for the analytic ansatz and an unacceptably 
large value of 5(1) for the NLO ChPT with finite volume corrections. In the future, it will be 
very interesting to see how the different ansatze for the chiral extrapolation become constrained or 
invalidated as we perform simulations with even lighter masses. We point out that the difference 
in the results from the analyses using the finite- volume ChPT and analytic ansatze is much smaller 
for the other quantities studied in this paper than for f K . 
The main physical results of this study are: 

fn = 124(2)(5)MeV {Eq.flSU}; h = 149(2) (4) MeV {Eq.©}; 

=£ = 1.204(7) (25) {Eq.(|63])}; 

JK 

mp(2GeV) = (96.2 ±2.7) MeV {Eq. fl§5)}; mJf(2GeV) = (3.59 ±0.21) MeV {Eq. flM])}; 

[L m (2GeV)] 1 / 3 = 256(6) MeV {Eq.®}; 

r = 0.487(9)fm and n = 0.333(9) fin {Eq. ®}. (103) 

For convenience we also display the equation number where the results were presented earlier in 
this paper to help the reader find the corresponding discussion. All the results in Eq. (11031) were 
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obtained after reweighting the strange-quark mass to its physical value at each /3, and the renor- 
malized quark masses were obtained using non-perturbative renormalization with non-exceptional 
momenta as described in Section|VIJ The low-energy constants obtained by fitting our data to 
NLO chiral perturbation theory can be found in Sec. lVEl 

The configurations and results presented in this paper are being used in many of our current stud- 
ies in particle physics phenomenology, including the determination of the Bk parameter of neutral 
kaon mixing in the continuum limit 113411 . In parallel to these studies we are exploiting config- 
urations generated at almost physical pion masses on lattices with a large physical volume (~ 
4.5 fm) but at the expense of an increased lattice spacing. Preliminary results obtained for the 



meson spectrum and 

r 

presented in Refs. |48, 



decay constants and for AI = 3/2 K — > %% decay amplitudes were recently 



69Q . Having access to data with excellent chiral and flavor properties with a 



range of lattice spacings and quark masses makes this an exciting time indeed for studies in lattice 
phenomenology. 
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Appendix A: Separate fits to 24 3 and 32 3 data 

In this section we report on results obtained by fitting the data from the 24 3 runs at /3 = 2. 13 and 
from the 32 3 runs at j6 = 2.25 separately to the predictions of SU(2) x SU(2) ChPT. This comple- 
ments the material presented in Sections [III] and [IV] in which we presented the results for masses 
and decays constants at each set of quark masses but did not perform the chiral extrapolations 
and also that in Section [V] in which we performed simultaneous chiral and continuum fits to the 
data at both lattice spacings. Our main motivation for studying separate fits here is to be able 
to compare directly our results obtained with the new data to those in our previous publication 
For that reason in this appendix we will be using the same renormalization constant as in 
our previous publication, which differs from the one used in the global analysis presented in the 
main part of this paper, see the discussion in Sec. [in] and App. |B]for details. We use the same 
method of iterated fits as outlined in our earlier publication [1]; at each lattice spacing we iterate 
the combined fits of the meson masses and decay constants with m x < 0.01 to the SU(2)-ChPT 
formulae, using kaon SU(2) ChPT to fit the kaon mass and decay constants and the extrapolation 
in the ^-baryon mass until convergence. The pion, kaon, and Q. masses are used to fix the phys- 
ical bare quark masses m u( i, m s and the lattice scale 1 /a. Predictions for the remaining physical 
quantities are then obtained by extrapolation to these physical quark masses. For further details 
see til]. In the case of the 24 3 ensembles, the runs have been extended since the publication of 
til] (see Sec. HI] and especially Tab. Ufor details) so that a direct comparison of the results from 
the previous (smaller) data set with the new extended data set is possible. We quote results from 
fits with and without corrections due to finite-volume effects. When including the finite volume 
corrections, the terms described in Appendix C of jlj] are included in the SU(2) ChPT in the pion 
sector (both for the meson masses and decay constants). We also include the correction terms 



containing the chiral logarithm of the light quark masses in the kaon decay constant [[820 and note 
that up to NLO in the light-quark masses, no finite- volume corrections arise in the masses of the 
kaon and ^-baryon. Below we present the physical results in the infinite-volume limit, i.e. after 
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removing the corrections. Finally, we will perform a naive continuum extrapolation of the results 
obtained by the separate fits at the two lattice spacings, which can then be compared to results 
from the combined chiral-continuum extrapolations using the global fits described in Sec.|V] Note 
that in this appendix also for the combined chiral-continuum extrapolations we are going to quote 
results obtained using our previous definition of Z^. For that reason the results reported here differ 
slightly from those in the main part of this paper. 



1. SU(2)-ChPT fits to 24 3 data 

In Tab. lXXXIXI we summarize our results from the iterative fits to the masses and decay constants 
measured on the 24 3 ensembles (see Sec. |In]for details) and compare them to our earlier results ob- 
tained with lower statistics yj]. We have performed two kinds of fits: one including the Q-baryon 
masses determined at all the simulated light-quark masses, m/ = 0.005, 0.01, 0.02, and 0.03, (as 
was done originally) and one where only the Q-baryon masses at the two lightest dynamical quark 
masses m/ = 0.005 and 0.01 are included. The latter, limited range is also the one used in the 
combined chiral-continuum extrapolations in Section |V] and in the separate fits to the 32 3 data 
in the next subsection. In Fig. [55] we plot the combined SU(2) ChPT fits (without finite-volume 
corrections) to the meson masses and decay constants in the pion sector. It is evident that over 
the fit range (m x + m y )/2 < 0.01, corresponding to a maximum meson mass of about 420 MeV, 
the data is well described by SU(2) ChPT. This is also true for the fits including the finite-volume 
corrections (not shown). 

We note that by comparing the results in the first two columns of Tab. IXXXIXI which have been 
obtained using the same (large) mass-range for the chiral extrapolation of the ^-baryon mass, the 
results obtained with the increased statistics (for each dynamical light-quark mass the statistics 
has nearly been doubled, see Section HITl) nicely agree with those from our previous publication yj] 
within the statistical uncertainty. Furthermore, we observe the expected reduction in the statistical 
error. For the remainder of the discussion, we focus on the fits in which only the two lightest 
dynamical masses have been included in the extrapolation of the ^-baryon mass, i.e. the last two 
columns of Tab. IXXXIXI The major difference resulting from this change in the fit range is in the 
value of the lattice scale I /a, but within 1.4 standard deviations (statistical error only, taking into 
account correlations) the results still show agreement. Including the finite-volume effects results 
in higher values for the decay constants (both in the chiral limit and at the physical point), which 
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Allton et al. [1J 
no FV-corr. 
£1: mi < 0.03 



increased statistics 
no FV-corr. incl. FV-corr. 

£1: mi < 0.03 £1: m t < 0.01 £1: m x < 0.01 



l/a [GeV] 

B m (2GeV) [GeV] 
/ [MeV] 

h 

f% [MeV] 
/* [MeV] 

fie/ fx 

mjjf (2GeV) [MeV] 
mp(2GeV) [MeV] 
m ud : 



1.729(28) 

2.52(0.11)(0.23) r 
114.8(4.1) 
3.13(0.33) 
4.43(0.14) 

124.1(3.6) 
149.6(3.6) 
1.205(0.018) 

3.72(0.16)(0.33) r 
107.3(4.4) (9.7) re 
1:28.8(0.4) 



1.731(19) 



1.784(44) 



1.784(44) 



2.63(0.06) (0.07) ren 2.69(0.09) (0.08) ren 2.63(0.09) (0.08) r 

111.5(2.9) 114.8(4.0) 117.1(4.0) 

2.76(0.24) 2.82(0.24) 2.59(0.27) 

4.54(0.10) 4.61(0.10) 4.57(0.11) 



121.2(2.5) 
147.9(2.6) 
1.220(0.011) 

3.56(0.08)(0.10) r 
101.0(1.9)(2.9) re 
1:28.37(0.27) 



124.4(3.6) 
151.0(3.7) 
1.214(0.012) 

3.48(0.12)(0.10) r 
99.0(3.0) (3.0) rei 
1:28.44(0.26) 



126.4(3.6) 
152.1(3.7) 
1.204(0.012) 

3.55(0.12) (0.1 l) r 
98.8(3.0)(3.0) rei 
1:27.89(0.28) 



aB 


2.414(61) 


2.348(43) 


2.349(44) 


2.298(45) 


af 


0.0665(21) 


0.0644(14) 


0.0643(14) 


0.0656(14) 


Lf x 10 4 


1.3(1.3) 


2.2(0.9) 


2.5(0.9) 


2.2(0.9) 


Lf x 10 4 


5.16(0.73) 


5.00(0.47) 


5.50(0.47) 


5.36(0.48) 


(24 2) -Lf)xl0 4 


-0.71(0.62) 


-0.09(0.45) 


0.03(0.45) 


0.01 (.49) 


(24 2) -4 2) ) x 104 


4.64(0.43) 


4.86(0.30) 


4.36(0.38) 


5.34(0.33) 


arhud 


0.001300(58) 


0.001331(43) 


0.001251(71) 


0.001274(72) 


am s 


0.0375(16) 


0.0377(11) 


0.0356(19) 


0.0355(19) 



TABLE XXXIX: Results from the SU(2) ChPT fits to the 24 3 data (without and with finite-volume correc- 
tions) compared to those from ||l|] obtained with lower statistics (without finite-volume corrections). We also 
quote in the lower part of the table the SU(2) ChPT fit parameters aB, af, Lj 2 ' (at the scale = 1 GeV) 
and bare quark masses am ul i s in lattice units. Only statistical uncertainties are quoted except for quark 
masses and the LEC B renormalized in the MS -scheme at 2 GeV where also the systematic uncertainty 
from the renormalization constant is quoted. (Mass renormalization constant at l/a = 1.73 1(19) GeV: 
Z m = 1.546(0.002) stat (0.044) ren and at l/a = 1.784(44)GeV: Z m = 1.559(0.003) stat (0.047) ren .) 
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0.080 




0.110 



0.105 



m, = 0.01 





FIG. 55: Combined SU(2) ChPT fits (without finite-volume corrections) for the meson decay constants (left 
column) and masses (right column) on the 24 3 data set at mi = 0.005 (top row) and 0.01 (bottom row). Only 
points marked with circles, corresponding to the range (m x + m y ) /2 < 0.01 are included in the fits. 

is a statistically significant effects (taking the correlations into account). In Tab. IXLI we compare 
the decay constants and their ratio obtained from the separate fits with the corresponding results 
from the global analysis at the simulated, finite value of the lattice spacing (i.e. not extrapolated 
to the continuum, see Sec.|V]and especially Tabs. IXXXfl IXXXTTl IXXXiTTl but note the difference 
due to the use of our previous definition of here). We are reassured by the observed agreement 
between the results obtained using the global fits with those obtained using our previous strategy 
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fn [MeV] /^[MeV] 


no FV-corr. 24 3 , j8 =2.13 separate 

global 


124.4(3.6) 151.0(3.7) 1.214(0.012) 
123(2) 150(2) 1.215(0.009) 


32 3 , j8 =2.25 separate 
global 


120.4(1.9) 147.1(2.0) 1.222(0.007) 
121(2) 147(2) 1.222(0.006) 


incl. FV-corr. 24 3 , j8 =2.13 separate 

global 


126.4(3.6) 152.1(3.7) 1.204(0.012) 
126(2) 151(2) 1.204(0.009) 


32 3 , j8 = 2.25 separate 
global 


122.3(1.9) 148.1(2.0) 1.212(0.007) 
123(2) 149(2) 1.210(0.006) 



TABLE XL: Comparison of the pion and kaon decay constants and their ratios at finite lattice spacing from 
separate (see Tabs. IXXXIXI IXLlT) and global fits using our previous definition of Z^. 

in Ref. [1] which was developed at that time to describe data at only a single lattice spacing. 



2. SU(2)-ChPT fits to 32 3 data 

The results of a separate fit on the 32 3 data set are summarized in Tab. IXLIl Here we only included 
the Q-baryon masses from the mi = 0.004, 0.006, and 0.008 ensembles. In Fig.[56]we show the fits 
for the meson masses and decay constants in the pion sector (without finite-volume corrections). 
Again, over the fit range ((m x + m y )/2 < 0.008), corresponding to a maximum pion mass of about 
400 MeV, the data is well described by SU(2) ChPT. 

As was already the case for the 24 3 ensembles, taking finite- volume corrections into account also 
leads to a good description of the data and results in higher values for the decay constants at the 
physical point and in the chiral limit. Again, taking the correlations into account, we note that this 
is a statistically significant effect. As was also the case on the 24 3 ensembles, we observe a good 
agreement for the decay constants and their ratio between the results of the separate fits to the 32 3 
data and the results from the global fits at finite lattice spacing, see Tab. IXLl 
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no FV-corr. 



FV-corr. incl. 



I /a [GeV] 

B m (2GeV) [GeV] 
/ [MeV] 
h 

u 

fx [MeV] 
f K [MeV] 
hlfn 

m^(2GeV) [MeV] 
miP(2GeV) [MeV] 
m u d : m s 



2.221(29) 



2.221(29) 



2.62(0.05) (0.06) ren 2.57(0.05) (0.06) r 
111.4(2.2) 113.7(2.2) 



2.84(0.21) 
4.18(0.09) 

120.4(1.9) 
147.1(2.0) 
1.222(0.007) 



2.61(0.24) 
4.10(0.09) 

122.3(1.9) 
148.1(2.0) 
1.212(0.007) 



3.58(0.07) (0.08) ren 3.64(0.07) (0.08) re 
100.6(1.7)(2.2) ren 100.4(1. 7) (2.2) rer 



1:28.08(0.19) 



1:27.60(0.20) 



aB 


1.826(0.024) 


1.790(0.025) 


af 


0.0502(0.0007) 


0.0512(0.0007) 


Lf x 10 4 


-0.75(0.79) 


-1.21082) 


Lf x 10 4 


5.14(0.40) 


4.87(0.41) 


(24 2) -Lf)xl0 4 


-0.93(0.42) 


-1.03(0.45) 


(24 2) -4 2) ) x 10 4 


6.22(0.23) 


7.37(0.24) 


arhud 


0.001040(31) 


0.001057(32) 


am s 


0.0292(08) 


0.0292(08) 



TABLE XLI: Results from the SU(2) ChPT fits to the 32 3 data (without and with finite-volume corrections). 
We also quote in the lower part of the table the SU(2) ChPT fit parameters aB, af, L) (at the scale 
K x = 1 GeV) and quark masses am ul i,s in lattice units. Only statistical uncertainties are quoted except for 
quark masses and the LEC B renormalized in the MS-scheme at 2 GeV where also the systematic uncertainty 
from the renormalization constant is quoted. (Mass renormalization constant at l/a = 2.221(29)GeV: 
Z m = 1.550(0.002) stat (0.034) ren .) 
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FIG. 56: Combined SU(2) ChPT fits (without finite-volume corrections) for the meson decay constants (left 
column) and masses (right column) on the 32 3 data set at mi = 0.004 (top row), 0.006 (middle row), and 
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3. Extrapolation to the Continuum Limit 

With the results obtained from separate chiral extrapolations on the 24 3 (extended statistics) and 
the 32 3 data sets (see the two previous subsections, respectively) we can perform a naive contin- 
uum limit extrapolation assuming a 2 -scaling. Of course, with only two lattice spacings available, 
we are not able to confirm this scaling behaviour. Further caveats include the fact that here, for 
simplicity, we did not use reweighting and so the dynamical strange-quark mass is not tuned to 
exactly the same value on the two data sets and indeed is not exactly the physical one on either 
set. Also, the dynamical light-quark mass ranges are a little different at the two lattice spac- 
ings, corresponding to unitary pion masses in the range 330-420 MeV on the coarser 24 3 lattices 
and 290-400 MeV on the finer 32 3 lattices (a similar statement is true for the partially-quenched 
masses). One might therefore expect a larger uncertainty in the chiral extrapolation of the 24 3 
results. In the naive continuum ansatz followed here, we are not taking into account this effect. 
Because of this, and maybe more importantly, since two separate chiral extrapolations have been 
performed (one at each of the two values of the lattice spacing), the continuum extrapolation is 
not completely disentangled from the chiral extrapolation. Recall that in our procedure for the 
global fits described in the main part of this paper, these two extrapolations are indeed disentan- 
gled. There this is achieved by adding 0(a 2 ) terms into the two functions, such that the chiral and 
continuum extrapolations are performed simultaneously and independently from each other. 
In Tab. IXLIII we repeat the results obtained at the two different lattice spacings (with and with- 
out finite-volume corrections) and give the values extrapolated to the continuum limit assuming 
a 2 scaling. Fig. 1571 illustrates the continuum extrapolation of the various quantities (only results 
obtained without taking into account finite-volume corrections are shown there). Note, that the 
two points at the different lattice spacings are completely uncorrelated, the only correlation in the 
data for the continuum extrapolation is between the uncertainty in the lattice spacing (the "x"- 
datum) and the quantity itself at that lattice spacing (the "y"-datum). These correlations were 
treated by the super-jackknife method which we have been using in our earlier work and which is 
clearly explained in [73, 74}] . For comparison, Tab. IXLlJ also contains our results from the com- 
bined continuum-chiral extrapolation as described in the main part of this paper but here using our 
previous definition of As one can see, the combined continuum-chiral extrapolation gives a 
substantially smaller (up to a factor of 5) statistical uncertainty compared to the naive continuum 
extrapolation. The main reason, of course, is the correlation in the combined fits between the two 
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no FV-corr. 






separate fits 


naive CL 


comb. chiral/CL 




24 3 ,/3 =2.13 32 3 , J8 =2.25 






a [fm] 


All C\f.t r )H\ 
V. 1 lUD^Z / ) 


A AOQQ/1 0\ 
V.Vooo(lZ) 


^0 


^0 


/ [MeV] 


114.8(4.0) 


111.4(2.2) 


105.2(10.4) 


107(2) 


h 


2.82(0.24) 


2.84(0.21) 


2.87(0.74) 


2.81(0.16) 


u 


4.61(0.10) 


4.18(0.09) 


3.39(0.36) 


3.76(0.08) 


fn [MeV] 


124.4(3.6) 


120.4(1.9) 


113.0(9.5) 


117(2) 


h [MeV] 


151.0(3.7) 


147.1(2.0) 


139.9(9.6) 


144(2) 


ft /fx 


1.214(0.012) 


1.222(0.007) 


1.236(0.030) 


1.233(0.008) 






includin 


g FV-corr. 






separate fits 


naive CL 


comb. chiral/CL 




24 3 ,/3 = 2.13 32 3 , J8 =2.25 






a [fm] 


0.1106(27) 


0.0888(12) 


^0 


^0 


/ [MeV] 


117.1(4.0) 


113.7(2.2) 


107.4(10.3) 


110(2) 


h 


2.59(0.27) 


2.61(0.24) 


2.64(0.83) 


2.55(0.18) 


u 


4.57(0.11) 


4.10(0.09) 


3.26(0.38) 


3.83(0.09) 


fn [MeV] 


126.4(3.6) 


122.3(1.9) 


114.8(9.4) 


119(2) 


h [MeV] 


152.1(3.7) 


148.1(2.0) 


140.9(9.6) 


145(2) 


hi fn 


1.204(0.012) 


1.212(0.007) 


1.226(0.029) 


1.219(0.007) 



TABLE XLII: Selected results from separate fits to the 24 3 and 32 3 data sets (Q. masses from m/ < 0. 1 for 
24 3 data set, cf. Tabs. IXXXIXI and IXLIb and their naive continuum limit assuming a 2 -scaling (see Fig. [57]) 
compared to results from the combined chiral-continuum extrapolation using the previous definition of Z^. 
The top table contains results without finite- volume corrections whereas the results in the bottom table were 
obtained by including finite-volume effects. 
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f [MeV] 
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FIG. 57: Results from separate fits (without finite-volume corrections) to the 24 3 and 32 3 data sets (black 
points) and the naive continuum-limit extrapolation (blue asterisks) for selected quantities assuming en- 
sealing. For details see Subsec. lA~3l and Tab. IXLIll 
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data sets at different lattice spacings. This correlation occurs because we require the fitted param- 
eters to be the same on both data sets and only include 0(a 2 ) corrections for the leading-order 
terms, as is consistent with our power counting scheme. In this way, the continuum extrapo- 
lation in the combined fits is also more constrained, leading to a smaller statistical uncertainty. 
Comparing the results of the naive continuum extrapolation and the combined continuum-chiral 
extrapolation for the quantities in Tab. IXLlJ we observe agreement better than 0.5-c (taking into 
account correlations) for all quantities except for I4, where the agreement still holds at the 1- or 
1.5-a level (without and with taking FV-corrections into account, respectively). It is reassuring, 
that the results from the two methods agree well, although the value of this statement is limited, 
given the large (statistical) uncertainty of almost 10% for the decay constants or even more in 
case of the LECs from the naive method. However, it should be noted that the same agreement 
holds, not only for the continuum values, but also for the results obtained in the separate fits as 
compared to the predictions of the global fit made for the finite lattice spacings. This has already 
been discussed in the previous subsections and is shown in Tab. IXLI 



Appendix B: Determining Z & 



As pointed out by Sharpe QJ/ZD and refined in Ref. [1], the normalization of the partially conserved 



axial current defined for domain wall fermions [75] is expected to deviate from that of the con 



ventionally normalized continuum current by an amount of order m KS a. Here and below when 
making such estimates we will introduce the explicit lattice spacing a and express the residual 
mass in physical units in order to make the comparison of various terms in a Symanzik expan- 
sion in powers of a easier to recognize. Since such a deviation can be viewed as 0(ma) which 
is formally larger than the 0(ma 2 ) which we neglect in our power counting scheme and because 
the normalization of this axial current plays a central role in our determination of the important 
quantities f n and fx, we have calculated this normalization factor Z^ numerically. We explain our 
method and result in this appendix. The first subsection contains a discussion of the theoretical 
issues and explains the basis for our method of determining Z^. The second subsection describes 
the actual calculation and results. 
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1. Determining the normalization of ^ 

To determine the normalization of ^ we compare the matrix element of four distinct domain wall 
fermion currents. The first two are the conserved/partially conserved vector and axial currents 
y^{x) and stf^{x) respectively, where a and /i are flavor and space-time indices. These currents 



were introduced by Furman and Shamir 117511 and involve fermion fields evaluated on each of the L s 
4-dimensional hyperplanes and at both the space-time points x and x + e^ where is a unit vector 
pointing the \l th direction. Thus, these currents are local but distributed in the fifth dimension and 
one-link non-local in space-time. While this vector current is exactly conserved, the divergence of 
the axial current contains the usual mass term and a mid-point term 7| ? . In the long-distance limit 
this midpoint term can be decomposed into the residual mass term, a piece that is conveniently 
written as (1 — Zg/) times the divergence of the same axial current and a final term of dimension 
five which we write out explicitly as the sum of the dimension-five, chiral rotation of the usual 
clover term and the four-dimensional Laplacian applied to the pseudoscalar density: 

Jl q = m KS q^X a q + ^^A^ + Cl qO^F^rq + C2d^q^X\ (Bl) 

In Equation (IB 1 1) X a is the generator which acts on the fermion fields corresponding to the flavor 
index a while q(x) and q(x) are the "physical", four-dimensional quark fields obtained by evaluat- 
ing the five-dimensional domain wall fields on the s = and s = L s — 1 boundaries. (See Eqs. (11) 
and (12) in Ref. QD.) 

The second pair of currents which we will need in this appendix is the local vector and axial 
currents, V£(x) and A^(jc), constructed in the standard way from the four-dimensional quark fields, 
q(x) and q(x). These currents are localized in all five dimensions and neither is conserved. 
Finally it will also be convenient to introduce the scalar densities q(x)q(x), q(x)X"q(x) from which 
the domain fermion mass is constructed and their chiral transforms q(x)^q(x), q(x)X a 'f'q(x). 
These four classes of operators will be labeled S(x), S a (x), P(x) and P a {x). 
Following Symanzik, we can add improvement terms to each of these six operators to insure that 
their Green's functions, when evaluated with an appropriately improved action, will agree with 
the corresponding continuum Green's functions up to errors of order a n . For our present purposes, 
accuracy up to 0(am) where m is a quark mass in physical units, will be sufficient. Since m res 
and m have a similar size, we are explicitly attempting to control the m res a corrections described 
above. We do not attempt to explicitly remove 0(a 2 ) terms since these will be eliminated by the 
final linear extrapolation a 2 — > 0. 
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In the discussion to follow we will recognize constraints on the required Symanzik improvement 
terms and relations between the various renormalization constants by applying the approximate 
chiral symmetry of domain wall fermions to Green's functions containing these various operators. 
For such arguments to be valid we will assume that these Green's functions are evaluated at suf- 
ficiently small distances that the effects of the vacuum chiral symmetry breaking of QCD can be 
ignored but at sufficiently large distances that the Symanzik improvement program can be applied. 
Since this discussion is a theoretical one, constraining the form of the Symanzik improvement 
terms, we need not be concerned about practical questions regarding the degree to which such 
conditions can be realized in our present calculation. 

Using the notation V^ a , A^ fl , S Sa and P Sa for the Symanzik-improved vector current, axial cur- 
rent, scalar density and pseudoscalar density respectively, keeping improvement terms which are 
nominally of order a and imposing charge conjugation symmetry, we find: 



y; fl = z r y^+c r d v qa^ v x a q (B2) 

Al a = Z^ + C^P" (B3) 

V* a = Z v V^ + C v d v qa^ y X a q (B4) 

A% = Z A A a ^+C A d^P a (B5) 

S Sa = Z s S a (B6) 

P Sa = Z P P a . (B7) 



In contrast to the Symanzik-improved current operators, we have not specified a normalization 
convention for the operators S Sa and P Sa . Adopting definitive conventions for S Sa and P Sa is not 
needed here beyond the requirement that those conventions are consistent with S Sa ±P Sa belonging 
to the (3, 3) / (3, 3) representations of the SU(3)lxSU(3)^ flavor symmetry. 
Because the operators S and P contain no vector indices, any correction terms must increase the 
dimension by two and we have chosen to neglect such 0(a 2 ) contributions. Thus, Eqs. (|B6I) and 
(IB7I) are particularly simple. However, we can also drop the dimension four, 0(a) correction terms 
to Eqs. (IB2I) - (IB5I) . This can be established by considering the chiral structure of the Symanzik 
and conserved/partially conserved current operators. Ignoring effects of order m, the Symanzik 
currents will couple to pairs of quarks which are either left- or right-handed. Likewise the domain 
wall conserved/partially conserved current operators couple to a pair of quarks with the same 
value of the coordinate s in the fifth dimension. For s = these are left-handed fermions while for 
s = L s — 1 they are right-handed. As the coordinate s moves into the fifth-dimensional bulk, the 
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amplitude for coupling to such physical modes decreases until when s ~ L s /2 the amplitude will 
be suppressed by two traversals half-way through the fifth dimension which implies a suppression 
of order m res a. Of course, the s ~ and s ~ L s — 1 terms will dominate. The character of the 
local vector and axial currents is simpler since they contain quark field strictly limited to s = 
and L s — 1. Since the four, dimension-four improvement terms included in Eqs. (1B2J) - (1B5I) involve 
pairs of quarks with opposite handedness, such terms require a complete propagation across the 
fifth dimension if they are to couple to the conserved/partially conserved or local currents. This is 
true even for the terms with general s which appear in the former currents. Thus, these correction 
terms involve an additional power of m KS a and are of order m res a 2 and can be neglected in our 
power counting scheme. 

With this simplification, we can demonstrate that to this order the following relations hold: 

Z r = 1 (B8) 
Zv = Z A (B9) 
Z s = Z P . (BIO) 

Equation (IB 81) follows easily from the fact that V£ is conserved at finite lattice spacing and has 
been given the conventional normalization. Equations (lB9b and (IB 101) can each be shown using 
essentially the same argument which we will now review. 

In the massless continuum theory the operators q c X a y^{ \ ±y 5 )q c are independent involving only 
right-handed/left-handed degrees of freedom. Here the label c indicates continuum. This implies 
the vanishing of the Symanzik-improved Green's function: 

(V^+Af)(x)(V^-A s v a )(y))=0. (Bll) 

This same property is obeyed by the local domain wall currents up to order (m KS a) 2 since non- 
vanishing terms which can contribute to the DWF version of Eq. (IB lib must connect both fermion 
degrees of freedom between the left and right walls requiring two-traversals of the fifth dimension 
and hence are of order (m res a) 2 [17,Q]. It is then easy to see that these two behaviors can be 



consistent through order m res a only if Zy = Za through order m ies a. We need only examine the 
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mixing between V„ a ±A^ a that is generated by Zy — Z A : 



((v^+A s p(x)-(v^-A s v a )(y)) 

= ((Z K y M fl +Z A A^)W-(Z v \/ v fl -Z A A fl v )(3;)) 

= \([{zy +z A )(v M fl +a«)(x) + (z y -z A )(y; -a; ) (*) 



(B12) 



• ' (Zy + Z A ) (V« -A a v ) (y) + (Z V - Z A ) (V v fl + A%) (y) ] ) 



The product of the left-most operators in the square brackets on the right-hand side of Eq. (IB 121) 
cannot mix at order m KS because of their construction from domain wall quark fields as explained 
above. Likewise the product of the right-most terms also vanishes. However, the two cross terms 
have non-zero correlators implying that for the entire expression to be of order m 2 es , the difference 
Zy — Za must be of order (m res a) 2 , demonstrating the intended result. A very similar argument 
can be constructed which shows that Z$ = Zp through order m ies a. One must invoke the flavor 
structure and, for example, consider correlators between (S 1 — iS 2 )(x) + (P 1 — iP 2 )(x)) and (S l + 
iS 2 )(y) + (P l + iP 2 )(y)) which also must vanish in the chiral limit. Here a = 1,2 is a specific 
choice of the eight octet indices a = 1 — 8. 

The relations in Eqs. (IB8I) . (IB9I) and (IB 101) were established by considering the domain wall and 
continuum theories in a limit in which the physical quark masses could be neglected, at sufficiently 
short distances that vacuum chiral symmetry breaking could be ignored but at sufficiently long 
distances that the Symanzik effective theory could be applied. While this is an excellent regime 
in which to establish these theoretical constraints, it is not a practical one for calculations. Thus, 
we will now employ these relations at low energies where vacuum chiral symmetry breaking is 
important in order to provide a practical method to compute Zjg. 

Since at low energies the left- and right-hand sides of Eqs. (|B4|) and (IB5I) must have identical 
matrix elements, the ratio of long-distance correlators computed with the Symanzik and local 
currents must give identical constants: Zy =Za- Thus, we have established: 



where we have introduced the fixed spatial index i, the temporal index and sources Vfiy) and 
P a (y) that will correspond to those used in our actual calculation. Next we can use the long- 



(Vf'Wiy)) (A s Q a (x)P a (y)) 
(V?(x)Vf(y)) " " (A a (x)P"(y)) 



(B13) 



119 



distance equality represented by Eqs. (IB2I) and (IB3I) to write 

(Ajj fl (jc)P fl (y)) 

Z ^ " «(^M>- (B15) 

Then we can combine Eqs. (IB 13b . (IB 141) and (IB 151) to yield an equation for Z^ which does not 
involve the Symanzik currents: 

_ (Ag(*)/*(y)) QTOjgbO) 
* «(x)J*(y)) <V?(*)V; a (y)> ' 

which determines Z^ in terms of four correlators which we have evaluated directly in our lattice 
calculation. 

In order to relate the discussion of the Symanzik improved operators given in Eqs. (IB2|) - (IB7I) with 
the operators appearing in Eq. (IB lb . we should recognize that the quantity Z^j has been introduced 
in two places. The most important is in the relation between the Symanzik current and the partially 
conserved domain wall operator in Eq. (IB3I) . It is this quantity that is determined in Eq. (IB 161) and 
which is needed to give a physical normalization to the axial current matrix elements determined in 
our calculation. However, the quantity Z^ also appears in the expression for J$ q given in Eq. (IB 1 b . 
For completeness, we will now demonstrate that these two quantities are in fact the same up to 
order (m Tes a) 2 . 

This is easily done by introducing a flavor-breaking mass term qMq into the DWF action, exam- 
ining the divergence equations obeyed by ¥2 and g/£ and using the relation Z$ = Zp established 
above. With the additional mass term the conserved/partially conserved vector and axial currents 
obey the lattice divergence equations, through 0(m ies a): 

tPV" = q[X a ,M]q (B17) 
tf^a = q{X a ,M}^q + 2m ies q^q-(Z^-\)A^^. (B18) 

Taking the Z^ — 1 term to the left hand side and recognizing that the scalar and pseudoscalar 
operators S a and P a are symmetrically normalized (Z$ = Zp), we can conclude that the operators 
Yn and Z^stf^ must be related to the corresponding Symanzik currents by the same factor. This 
establishes that our two definitions of Z^ are consistent. 

We will conclude this analysis with a brief discussion of the effects of the explicit quark mass, 
rrif, on the operator product expansion represented by Eq. (IB 1 b and on the Symanzik-improved 
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operators given in Eqs. (IB2I) - (IB7I) . Although nif explicitly connects the s = and s = L s — 1 
walls, it can combine with the midpoint operator J^ q appearing on the left hand side of Eq. (IB II) 
to create effects with arbitrary chiral properties. Thus, we expect multiplicative corrections of 
the form (1 + £,m/a) i<,<4 to each of the four terms on the right hand side of Eq. (IB1|) . In the 
case of the left-most term the correction is of order mfin KS a while for the remaining three terms 
the corrections are of order m^m res a 2 or m^m res a 3 , all beyond the level of accuracy of the current 
paper. The conclusion that Zy = 1 through order m res a 2 (and order mya 2 ) prevents the appearance 
of a factor 1 +b(mfa) multiplying the Z-y in Eq. (IB2I) . The argument that Z& = Z v and Z$ = Sp 
with corrections of order (m Tes a) 2 applies equally well to the left-right mixings created by nif but 
again the allowed mftn ies a 2 and (nif a) 2 terms are negligible within our present power counting 
scheme so Eqs. (|B4|) - (IB7I) need no 0(mfa) corrections. Lastly, consider adding a factor of the 
form (1 +b(mfa)) multiplying the Z^ on the right-hand side of Eq. IB3I As explained above, a 
similar correction to Z^ appearing in Eq. (IB 1 1) carries the additional suppression of one power 
of m KS a. Since the equality derived above between the Z^ factors appearing in the divergence 
equation, Eq. (IB1|) . and the Symanzik-improved current in Eq. (IB3I) . holds at order nif a such 
a 1 + b(mfa) factor is not allowed in Eq. (IB3I) . Thus, no nifa terms need to be introduced into the 
equations presented in this appendix. 



2. Computational method and results 

We have evaluated the two factors in Eq. (IB 161) to determine Z^ on both the 32 3 x 64, /3 = 2.25 
(mi = 0.004, 0.006 and 0.008) and the 24 3 x 64, /3 = 2. 13 (m/ = 0.005, 0.01 and 0.02) ensembles. 
We used a small subset of these six ensembles and obtained the results given in Tab. IXLIIII The 
results presented for Z^jZ^ duplicate those from the calculation of described in Sections III 
and IV. In this appendix we add the factor Z^ in the denominator because we are now determining 
the deviation of this factor from unity. We do not simply use the results presented earlier in the 
paper because our calculation of Zy /Z-y has been performed on a subset of the configurations 
analyzed earlier and results for Zj^jZ^ are needed on this same subset of configurations if ratios 
with meaningful jackknife errors are to be determined. 

The ratio Za/Z^ was computed from the same ratio of current-pseudoscalar correlators studied in 



Sections III and IV, using the method specified in Ref. llwh . Similar methods are used to compute 
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j8 mi Z A /Z^ Zy/Zy Zy/Z^ Fit range iV meas 

2.13 0.02 0.71900(20) 0.6956(17) 1.0336(25) 9-54/9-17 50 

2.13 0.01 0.71759(16) 0.6998(20) 1.0254(29) 9-54/9-17 50 

2.13 0.005 0.71743(30) 0.6991(17) 1.0262(25) 9-54/10-19 105 
2.13 -m res 0.71615(36) 0.7019(26) 1.0208(40) 

2.25 0.008 0.74526(12) 0.73802(55) 1.0098(7) 9-54/9-20 85 

2.25 0.006 0.74523(12) 0.73853(64) 1.0090(9) 9-54/9-18 76 

2.25 0.004 0.74513(15) 0.73871(77) 1.0087(10) 9-54/10-19 166 
2.25 -m res 0.74499(34) 0.7396(17) 1.0073(23) 

TABLE XLIII: Results for the ratios Z^/Z^, Zy /Zy and Zy jZ^ computed on six ensembles. The rows 
with quark mass — m res contain the chiral extrapolation to the light quark mass m\ = — m res . The left-hand 
portion of the fit range gives that used for the axial current ratio while the right hand portion that for the 
vector current. For the Zy /Zy calculation the data at t and 63 — t were combined for < t < 32. 

Zy /Zy using the ratio of vector correlators 

/ = - 7 ^— f, (B19) 

an equation expected to be valid for time separations / much larger than one lattice spacing: t^> a. 
Figure[58] shows the right-hand side of Eq. (IB 191) as a function of time for the case of the lightest 
mass for each of the 32 3 and 24 3 ensembles. A constant fit to plateau regions identified by the 
horizontal lines was then used to determine the Zy/Zy on the left-hand side of this equation. 
Fig. [59] displays the chiral extrapolation of the two quantities Z^/Z^ and Zy/Zy on both sets of 
ensembles. 

Two useful results follow from this Appendix. First the ratio Zy /Z^ differs from unity on our 
two ensembles and that difference decreases more rapidly than a 2 with increasing /3. Thus, we 
will obtain more accurate results in our continuum extrapolation from both matrix elements of 
the local axial current and our NPR calculations which are normalized using off-shell Green's 
functions containing the local vector and axial currents if we convert the normalization of these 
local currents to the usual continuum normalization by using the ratio Zy/Zy instead of the ratio 
Za/Z^, the quantity which we have used in previous work for such conversions. The values of 
Zy/Zy presented in Table IXLiTTl are therefore used to normalize the results presented in the current 
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FIG. 58: Plots of the correlator ratio which determines the renormalization factor Zy/Zy as a function of 
time. The left panel shows results from the 32 3 , mi = 0.004 ensemble while the right panel the result from 
the 24 3 , m\ = 0.005 ensemble. The horizontal line with error bands in each panel shows the fitting range 
and the result obtained in each case. 




am q am q 

FIG. 59: The quantities Za/Z^ and Zy jZy extrapolated to the chiral limit for the 32 3 (left panel) and 24 3 
(right panel) ensembles. 

paper and are the second result obtained in this appendix. Because these ratios were calculated 
on a smaller subset of configurations than were used for our main results, we have included their 
statistical fluctuations as independent within our superjackknife, statistical error analysis. Since 
these fluctuations are at or below the 0.5% level, this omission of possible statistical correlations 
is unimportant. 
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Appendix C: Statistical errors of re weigh ted quantities 

In this appendix we discuss the statistical errors that should be expected when Monte Carlo data 
is reweighted to obtain results for a gauge or fermion action that is different from that used to 
generate the data. Throughout this discussion we will make the assumption that the reweighting 
factors are not correlated with the data. Of course, if this assumption were exactly true then the 
reweighting would not be needed. However, the correlation between the data and reweighting 
factors is often small in practice and neglecting this correlation may well provide a reasonably 
accurate view of the resulting errors. As we will show, with this assumption the usual analysis of 
the statistical errors applies easily to reweighted data and yields simple, useful formula which we 
present here. 

Consider a quantity x and the corresponding ordered ensemble of N Monte Carlo configurations 
with corresponding measured values {x n }, 1 < n < N. For each of these N configurations we will 
determine a reweighting factor w n so that the final, reweighted quantity of interest is given by 

<*>* = %^. (CD 

Here the single brackets (. . .) N indicate an average over a single Monte Carlo ensemble of N 
samples. In this appendix we are interested in how the statistical fluctuations in the quantity (x) N 
are affected by the operation of reweighting. We can then express the true value for xn as 

xn = (((x) n )) (C2) 

where the double brackets ((...)) indicate a "meta" average over many equivalent Monte Carlo 
ensembles. The statistical fluctuation present in a particular result (x) N can then be characterized 
by the average fluctuation of (x) N about xj/: 



Error(;c) = 




(C3) 



A quantity such as (xn), defined in Eq. (IC1I) as a ratio of averages, will be a biased estimator 
of the physical result which must be determined in the limit N — > °°. Thus, the meta average 
xn = (((x)n)) Wli l differ from the true result by terms of order 1 /N. While these 1 /N corrections 
are not difficult to enumerate and estimate from our data, these corrections are not the subject of 
the present appendix and will not be considered further here. Instead we will study how the size 
of the statistical fluctuations of (x N ) about xjj is affected by the reweighting. Thus, the quantity 
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LnWn J y Ln'W n > 

(E„=l X n W n ~ X^Ln=l Wn) (£%=l *>W ~ */V Eff = l W^) 

(Ew=l(^n -^V)W»)) (E^ = i(-V -^v)>v) 




Error(jc) defined in Eq. (IC3I) describes the average deviation of from not from the iV — > °° 
limit of J/y. 

We will now work out an expression for Error(x) in the case that nearby measurements x„ and 
x n+ i in a single Markov chain (or reweighting factors w„ and w„+/) are correlated but with the 
assumption that x n and w n+ / are not: 

W,-^) 2 )) - // f%^-^ fe^W) )) (C4) 



(C5) 
(C6) 

(C7) 



where in the last line we have used our assumption of the lack of correlation between the x n and 
w n to write the average of their product as the product of their separate averages. We have also 
assumed that our sample size N is sufficiently large that correlated fluctuations of the averages in 
the numerator and denominator will be sufficiently small that the average of the original ratios and 
products can be replaced by the corresponding ratios and products of the individual averages. 
This result can be cast in a simple form if we define the three averages: 

5x 2 = (((x„-^) 2 )) (C8) 
w = ((w n )) (C9) 
w 1 = ((wl)\ (CIO) 



(where 8x 2 is the usual width of the distribution of the measured quantity x n ) and the two autocor- 
relation functions: 

' K X n — Xfj) (x n+ i — Xjj) 



a,) = * s? <C11) 



W n Wn+l 

W(l) = " _ " , (C12) 



w 2 



defined so that C(0) = W(0) = 1. Making the conventional assumption that the range of / over 
which the correlation function C(l) is non-zero is small compared to the sample size N and using 
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the quantities defined above, we can rewrite Eq. (IC7I) as 



2 \\ 8x 2 tfr\C(l)W(l)w 



2 



(x) N -m >> = — :; a ;_, 2 (ci3> 

/ // N(w) 
where the autocorrelation time T cor r is defined as 

^ma x 

Tcorr = £ C(l)W(l). (C15) 

* — ^max 

The limit L max is chosen to be larger than the region within which C(l) is non-zero and has been 
introduced as a reminder that when working with a single finite sample, one must take care to 
evaluate the limit of large N before the limit of large L max . Finally, Eq. (IC14I) can be written in the 
conventional form 

1 8 X 2 

ErrorW = V^ <C16) 

where the effective number of configurations N e ff is given by: 

N w 2 

Nc&=—- =■ (C17) 

Tcorr W 

This result makes precise a number of aspects of reweighting that may be useful to understand. 
In the case that there are no autocorrelations so T cor r = 1> the ratio w 2 /w 2 expresses the degree to 
which the reweighting process selectively samples the original data and degrades the initial statis- 
tics. The general inequality w 2 /w 2 < 1 (a consequence of the Schwartz inequality) is saturated 
only in the case that the reweighting factors w n do not vary with n. In the extreme case that a 
single sample w n dominates the averages then v? 2 /w 2 = l/N and A^ e ff = 1. Thus, in the case of 
uncorrected data (which is the case for most of the results presented here) we should expect the 
statistical fluctuations to grow as the degree of reweighting increases by the factor w 2 /w 2 . 
Including autocorrelations makes the effects of reweighting on the size of the statistical fluctu- 
ations less certain because the behavior of the factors l/T CO rr and w 2 /w 2 in Eq. (IC17I) become 
entangled. In the limit in which the autocorrelation time associated with the measured quantity x n 
alone, 



E C(0, (C18) 
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becomes much larger than that of the reweighting factor w n , then the majority of the sum in 
Eq. (IC 151) contributing to T corr will come from values of / where (lw n w n+ i^) « ((w^) so mat 

—2 

T corr « T x ^. (C19) 

In this case the error given by Eq. (IC16I) reduces to the standard expression a/ 8x 2 t x /N that holds 
if no reweighting is performed! Of course, this is easy to understand. When such long autocorre- 
lation times are involved, the average over the autocorrelation time is providing an average over 
the reweighting factors w n which is sufficiently precise that the error-enhancing fluctuations in 
the reweighting factors are averaged away. Given the large size of the fluctuations between the 
reweighting factors and the relatively short autocorrelation times seen in our data, it is unlikely 
that this averaging would be seen in the results presented here. 

A second type of behavior for T corr occurs if the w„ are relatively uncorrected and w 2 3> w 2 so that 
only the / = term contributes to the sum in Eq. (IC15I) giving T corr = 1 . In this case reweighting 
has removed the effects of autocorrelation but increased the statistical fluctuations by the factor 
w 2 /w 2 which was assumed to be large. Here the fluctuation-enhancing effects of autocorrelations 
and reweighting are not compounded. 
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